home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / C / Applications / UIFlow 1.0.1 / UIFlow Source / Kevin / uiflowvs.f < prev   
Encoding:
Text File  |  1992-04-23  |  176.9 KB  |  5,208 lines  |  [TEXT/????]
  1. C======================================================================C
  2. C                                                                      C
  3. C                      UIFLOW - 2 D                                    C
  4. C                    -----------------                                 C
  5. C                                                                      C
  6. C                                                                      C
  7. C     A computer program for TWO-dimensional analysis of               C
  8. C     turbulent elliptic flows in general coordinates.                 C
  9. C                                                                      C
  10. C     UIFLOW2D uses a decoupled multigrid procedure                    C
  11. C     and solves for cell-centered velocities.                         C
  12. C                                                                      C
  13. C     Changes made to this code:                                       C
  14. C     -------------------------                                        C
  15. C                                                                      C
  16. C     (1)  Inclusion of subroutine vset.                               C
  17. C     (2)  Specification of o2 mass fractions based on fuel fractions. C
  18. C     (3)  Error corrected in sub. gammod (-imaxl).                    C
  19. C     (4)  Adaptive cycle removed from sub. moment.                    C
  20. C     (5)  Correct placement of u and v residual calculations.         C
  21. C                                                                      C
  22. C======================================================================C
  23.       SUBROUTINE UIFLOW
  24.       include 'UIFlow.com'
  25. C
  26.       CEXTERNAL SHOWLINE
  27.       CEXTERNAL STOPKEY
  28.       CEXTERNAL ABORTC
  29.       INTEGER*2 STOPKEY
  30.       CHARACTER*255 buffer
  31.       CHARACTER*1 NULL
  32. C
  33.       open(4, file='UIFlow.plt' ,access='sequential',status='unknown')
  34.       open(5, file='UIFlow.In'  ,access='sequential',status='unknown')
  35.       open(8, file='UIFlow.out' ,access='sequential',status='unknown')
  36.       open(11,file='UIFlow.Grid',access='sequential',status='unknown')
  37. C
  38. C.... Read input data and initialise flowfields
  39. C
  40.       NULL = char(0)
  41.       call SHOWLINE('Beginning Computation')
  42.       call  param
  43.       call  header
  44.       call SHOWLINE('Reading Input File')
  45.       call  input
  46.       if ( model .eq. 2  ) call search
  47.       call  outp
  48.       call  const
  49.       if ( model .ne. 0  ) call constr
  50.       call  zero
  51.       if ( kpgrid .eq. 1 ) call gridp
  52.       if ( kpgrid .eq. 0 ) call grid
  53.       call  geomm
  54.       call  init
  55. C
  56.       write (8,1001) 
  57.       write (8,1002)
  58. C
  59. C.... Set multigrid cycle
  60. C
  61.       igrf = 1
  62.       wrkunt = 0.0
  63.       call SHOWLINE('Writing History')
  64.       write (8, 1006) 'Convergence History'
  65.       write (8, 1007)
  66.       write (8, 1008) igrf
  67.       write (8, 1013)
  68.       write (8, 1005)
  69. C
  70.       write (buffer, 2013) NULL
  71.       call SHOWLINE(%ref(buffer))
  72.       write (buffer, 2008) igrf, NULL
  73.       call SHOWLINE(%ref(buffer))
  74.       write (buffer, 2013) NULL
  75.       call SHOWLINE(%ref(buffer))
  76.       write (buffer, 2005) NULL
  77.       call SHOWLINE(%ref(buffer))
  78.       write (buffer, 3006) NULL
  79.       call SHOWLINE(%ref(buffer))
  80.       write (buffer, 3007) NULL
  81.       call SHOWLINE(%ref(buffer))
  82. C
  83.       xxx = 0
  84. 10    continue
  85. C
  86.       nitn(igrf) = nitn(igrf) + 1
  87.       if ( model .ne. 0 ) call lamvis(igrf)
  88.       if ( klam  .eq. 0 ) call tvis  (igrf)
  89. C
  90. C.... Solve momentum and continuity equations
  91. C
  92.       call  moment(igrf)
  93. C
  94. C.... Solve scalar equations including k and epsilon
  95. C
  96.       call  scalrs(igrf)
  97. C
  98. C.... Update thermodynamic properties.
  99. C
  100.       if ( (model .eq. 1) .or. (model .eq. 2) ) then
  101.        call props (igrf)
  102.       endif
  103. C
  104. C.... Enforce boundary conditions.
  105. C
  106.       call  extrpl(igrf)
  107. C
  108. C.... Check convergence on fine grid
  109. C
  110.       write(8,1009)  igrf, nitn(igrf), (error(igrf,nv), nv=3,8)
  111.       write(buffer,2009)  igrf, nitn(igrf), (error(igrf,nv),nv=3,8), NULL
  112.       call SHOWLINE(%ref(buffer))
  113.       xxx = STOPKEY()
  114.       if (xxx .eq. 1) goto 9999
  115. C
  116.       if ( error(igrf,3) .lt. tolr(igrf) ) goto 20
  117.       if ( nitn(igrf)    .ge. maxitn     ) goto 9000
  118.       if ( error(igrf,3) .gt. 1.0E+5     ) goto 9000
  119.       go to 10
  120. C
  121. 20    continue
  122.       write (8,1010) igrf
  123.       write (8,1012)
  124. C
  125.       write (buffer,2012) NULL
  126.       call SHOWLINE(%ref(buffer))
  127.       write (buffer,2010) igrf, NULL
  128.       call SHOWLINE(%ref(buffer))
  129.       write (buffer,2012) NULL
  130.       call SHOWLINE(%ref(buffer))
  131. C
  132. C.... Terminate on finegrid or prolongate
  133. C
  134.       if(igrf .eq. ngrid) go to 50    
  135. C
  136. C.... Extrapolate Solution and increment grid
  137. C
  138.       call  prlong(igrf)
  139.       igrf = igrf + 1
  140.       call  extrpl(igrf)
  141.       write (8,1006) 'Convergence History Continued'
  142.       write (8,1007)
  143.       write (8,1008) igrf
  144.       write (8,1013)
  145.       write (8,1005)
  146. C
  147.       write (buffer,2008) igrf, NULL
  148.       call SHOWLINE(%ref(buffer))
  149.       write (buffer,2013) NULL
  150.       call SHOWLINE(%ref(buffer))
  151.       write (buffer,2005) NULL
  152.       call SHOWLINE(%ref(buffer))
  153.       go to 10
  154. C
  155. 50    continue
  156. C
  157. C.... Converged on the desired finest grid
  158. C
  159.       write (8,1015) 'The Problem Has Converged !!'
  160.       write (8,1020)
  161. C
  162.       write (buffer,2020) NULL
  163.       call SHOWLINE(%ref(buffer))
  164.       write (buffer,2015) 'The Problem Has Converged !!  Wait For Final
  165.      1 Results to be Printed to Output File!', NULL
  166.       call SHOWLINE(%ref(buffer))
  167.       write (buffer,2020) NULL
  168.       call SHOWLINE(%ref(buffer))
  169.       go to 9001
  170. C
  171. 9000  continue
  172.       write (8,1015) ' The Program is Not Converging Well,
  173.      1 Check The Input Data'
  174. C
  175.       write (buffer,2015) ' The Program is Not Converging Well,
  176.      1 Check The Input Data', NULL
  177.       call SHOWLINE(%ref(buffer))
  178.       call ABORTC()
  179.       goto 9999
  180. C
  181. 9001  continue
  182.       call  extrpl(ngrid)
  183.       write (8,1006) 'Converged Flow Fields On Finest Grid'
  184.       write (8,1007)
  185. C
  186.       write (buffer,2007) NULL
  187.       call SHOWLINE(%ref(buffer))
  188.       write (buffer,2006) 'Converged Flow Fields On Finest Grid',NULL
  189.       call SHOWLINE(%ref(buffer))
  190.       write (buffer,2007) NULL
  191.       call SHOWLINE(%ref(buffer))
  192. C
  193.       call  output(ngrid)
  194.       call  vset
  195. C      call  ftout (ngrid)
  196. C
  197. 990   format (5x,'Reading Input Data')
  198. 1001  format (72('*')/15x,'Initial Fields For Coarsest Grid')
  199. 1002  format (72('*'))
  200. 1005  format (/3x,'Grid',2x,'Iter',3x,'Mass',5x,8('-'),
  201.      1  'Norm. Residuals in Scalars ',8('-')/3x,'Numb',
  202.      2  2x,'Numb',3x,'Residual',3x,'Swirl',4x,'Enthalpy',4x,
  203.      3  'k',4x,'eps',4x,'Mix. Fraction'/3x,'----',2x,'----',3x,
  204.      4  8('-'),3x,5('-'),4x,8('-'),2(3x,3('-')),4x,13('-'))
  205. 2005  format (3x,'Grid',2x,'Iter',3x,'Mass',5x,8('-'),
  206.      1  'Norm. Residuals in Scalars ',8('-'),a1)
  207. 3006  format (3x,'Numb',2x,'Numb',3x,'Residual',3x,'Swirl',4x,'Enthalpy',
  208.      1  4x,'k',4x,'eps',4x,'Mix. Fraction',a1)
  209. 3007  format (3x,'----',2x,'----',3x,8('-'),3x,5('-'),4x,8('-'),
  210.      1  2(3x,3('-')),4x,13('-'),a1)
  211. 1006  format (/72('*')/25x,72a)
  212. 2006  format (25x,72a,1a)
  213. 1007  format (72('*'))
  214. 2007  format (5x,72('*'),1a)
  215. 1008  format (72('-')/15x,'Starting Grid Number',i5)
  216. 2008  format (15x,'Starting Grid Number',i5,a1)
  217. 1009  format (5x,i2,3x,i4,3x,6(1pe12.5,1x))
  218. 2009  format (5x,i2,3x,i4,3x,6(1pe12.5,1x),a1)
  219. 1010  format (/5x,30('+')/6x,'Converged on Grid Number',i4)
  220. 2010  format (6x,'Converged on Grid Number',i4,a1)
  221. 1012  format (5x,30('+'))
  222. 2012  format (5x,30('+'),a1)
  223. 1013  format (72('-'))
  224. 2013  format (72('-'),a1)
  225. 1015  format (/5x,81('*')/5x,81a)
  226. 2015  format (5x,81a,a1)
  227. 1020  format (5x,81('*'))
  228. 2020  format (5x,81('*'),a1)
  229. 1028  format(i4,a1,1pe10.3)
  230. C
  231. 9999  continue
  232.       close(UNIT = 4, status='keep')
  233.       close(UNIT = 5, status='keep')
  234.       close(UNIT = 8, status='keep')
  235.       close(UNIT = 11,status='keep')
  236. C
  237.       return
  238.       end
  239. C======================================================================C
  240.        subroutine adjstx (igrl)
  241. C                                                                      C
  242. C      Purpose:  Perform integral mass adjustments to satisfy overall  C
  243. C                mass balance at every zi station.                     C
  244. C                                                                      C
  245. C======================================================================C
  246.       include 'UIFlow.com'
  247.       include 'UIFlow.indx'
  248. C
  249.       delp  = 0.0
  250.       sumin = 0.0
  251. C
  252.       do 10 j = 2, jmax1
  253.            ij = 1 + (j-1) * imaxl + ibeg
  254.        sumin  = sumin + cx(ij)
  255. 10    continue
  256. C
  257.       do 30 i = 2, imax2
  258.         ijf   = i + ibeg
  259.         ijl   = i + (jmax1-1) * imaxl + ibeg
  260.        sumin  = sumin + cy(ijf) - cy(ijl)
  261. C
  262. C.... Calculate actual mass flow rate.
  263. C
  264.        sumflx = 0.0
  265.        sumcf  = 0.0
  266. C
  267.        do 20 j = 2, jmax1
  268.             ij = i + (j-1) * imaxl + ibeg
  269.         sumflx = sumflx + cx(ij)
  270.         sumcf  = sumcf  + ae(ij)
  271. 20     continue
  272.        if (i .eq. imax1) sumcf = 1.0
  273. C
  274. C.... Scale fluxes to satisfy mass continuity, and change pressure to
  275. C     correspond to change in mass flux.
  276. C
  277.        delp    = delp + (sumin - sumflx) / sumcf
  278.        if (i .eq. imax1) delp = 0.0
  279.        do 25 j  = 2, jmax1
  280.             ij  = i + (j-1) * imaxl + ibeg
  281.         cx(ij)  = cx (ij) * sumin / sumflx
  282.         cu(ij)  = cu (ij) * sumin / sumflx
  283.         p(ij+1) = p(ij+1) - delp
  284. 25     continue
  285. 30    continue
  286. C
  287.       return
  288.       end
  289. C======================================================================C
  290.       subroutine apcalc (igrl)
  291. C                                                                      C
  292. C     Purpose:  Calculate apu and apv for use in computing             C
  293. C               the mass fluxes.                                       C
  294. C                                                                      C
  295. C======================================================================C
  296.       include 'UIFlow.com'
  297.       include 'UIFlow.indx'
  298. C
  299.       do 11 j   = 2, jmax1
  300.           ioff  = (j-1) * imaxl + ibeg
  301.        do 10 i  = 2, imax1
  302.              ij = i + ioff
  303.            ijp1 = ij + 1
  304.            ijm1 = ij - 1
  305.             ijp = ij + imaxl
  306.             ijm = ij - imaxl
  307.         apu(ij) = (ae(ij)  * u(ijp1) + aw(ij) * u(ijm1) +
  308.      1             an(ij)  * u(ijp)  + as(ij) * u(ijm)  +
  309.      2             apu(ij) ) / app(ij)
  310.         apv(ij) = (ae(ij)  * v(ijp1) + aw(ij) * v(ijm1) +
  311.      1             an(ij)  * v(ijp)  + as(ij) * v(ijm)  +
  312.      2             apv(ij) ) / ap(ij)
  313.         apm(ij) =  ap(ij)
  314. 10     continue
  315. 11    continue
  316. C
  317.       return
  318.       end
  319. C======================================================================C
  320.       subroutine bbndry
  321. C                                                                      C
  322. C     Purpose:  Prescribe boundary conditions on the eta minus         C
  323. C               (bottom) boundary of the calculation domain.           C
  324. C                                                                      C
  325. C======================================================================C
  326.       include 'UIFlow.com'
  327.       igrl = ngrid
  328.       include 'UIFlow.indx'
  329. C
  330.       j = 1
  331.       do 12 n = 1, nsym
  332. C
  333.        if (model .eq. 0) rhym(n) = rhgs
  334.        if (model .eq. 1) rhym(n) = pref * wmair / (gascon * tym(n))
  335.        if (model .eq. 2) then
  336.         yn2ym   = 1.0 - fuym(n) - o2ym(n) - h2oym(n) - co2ym(n)
  337.         rwmym   = fuym(n)/wmfu  + o2ym(n)/wmox + yn2ym/wmn2
  338.      1          + co2ym(n)/wmco2 + h2oym(n)/wmh2o
  339.         wmym(n) = 1.0 / rwmym
  340.         rhym(n) = pref * wmym(n) / (gascon * tym(n))
  341.        endif
  342. C
  343.           ifl = ifym (n,ngrid)
  344.       ill = ilym (n,ngrid)
  345.        do 10 i = ifl, ill
  346.          ijf   = i + (j-1) * imaxl + ibeg
  347.         u   (ijf)  = ubym(n)
  348.         v   (ijf)  = vbym(n)
  349.         w   (ijf)  = wbym(n)
  350.         p   (ijf)  = 0.0
  351.         t   (ijf)  = tym (n)
  352.         tke (ijf)  = tkym(n)
  353.         tde (ijf)  = tdym(n)
  354.         f   (ijf)  = fym(n)
  355.         g   (ijf)  = gym(n)
  356.         yfu (ijf)  = fuym(n)
  357.         yo2 (ijf)  = o2ym(n)
  358.         yh2o(ijf)  = h2oym(n)
  359.         yco2(ijf)  = co2ym(n)
  360.         yn2 (ijf)  = yn2ym
  361.         gam (ijf)  = vscty
  362.         amu (ijf)  = vscty
  363.         amut(ijf)  = cd * rhym(n) * tkym(n) * tkym(n) /
  364.      1               ( tdym(n) + 1.0e-30 )
  365.         wmol(ijf)  = wmym(n)
  366.         rho (ijf)  = rhym(n)
  367.         cv  (ijf)  = a21(ijf)*u(ijf) + a22(ijf)*v(ijf)
  368.         cy  (ijf)  = cv (ijf)*rho(ijf)
  369. 10     continue
  370. 12    continue
  371.       return
  372.       end
  373. C======================================================================C
  374.       subroutine bcor (igrl, q1)
  375. C                                                                      C
  376. C     Purpose:  Extrapolate interior values of the given quantity to   C
  377. C               the wall.  Consistent with enforcing a zero normal     C
  378. C               derivative of the given variable.                      C
  379. C                                                                      C
  380. C======================================================================C
  381.       include 'UIFlow.com'
  382.       dimension q1(*)
  383.       include 'UIFlow.indx'
  384. C
  385.       do 10 j = 1, jmaxl
  386.          ijf  = (j-1) * imaxl + ibeg
  387.        q1(ijf + 1)   = q1 (ijf  + 2)
  388.        q1(ijf+imaxl) = q1 (ijf+imax1)
  389. 10    continue
  390. C
  391.       do 20 i = 1, imaxl
  392.          ijf  = i + ibeg
  393.          ij1  = ijf + jmax1 * imaxl
  394.        q1(ijf) = q1 (ijf+imaxl)
  395.        q1(ij1) = q1 (ij1 - imaxl)
  396. 20    continue
  397. C
  398.       return
  399.       end
  400. C======================================================================C
  401.       block data param
  402. C======================================================================C
  403.       include 'UIFlow.com'
  404.       data cd, cdqtr, cdtqtr, ee, ak/ 0.09,0.5477,0.1643,9.025,0.4/
  405.       data ce1, ce2, cg1, cg2, cr/ 1.44, 1.92, 2.8, 2.0, 3.0/
  406.       data acpox, bcpox, ccpox/938.5681,0.0972,-1.7167e-5/
  407.       data acpn2, bcpn2, ccpn2/842.5130,0.2364,-6.1891e-5/
  408.       data acpco2,bcpco2,ccpco2/842.7291,0.2927,-7.8040e-5/
  409.       data acph2o,bcph2o,ccph2o/1215.6587,0.7182,-1.3888e-4/
  410.       data acpfu, bcpfu, ccpfu/1420.5564,1.7792,-3.9418e-4/
  411.       data wmh2o, wmco2/ 18.01520, 44.00980/
  412.       data wmn2, wmox/ 28.161 , 31.99880/
  413.       data hox , hn2 / 2.8802e5, 2.7057e5 /
  414.       data hco2, hh2o/ 2.7521e5, 4.2261e5 /
  415.       data rox, wmair, gascon/0.2331, 28.967, 8314.3/
  416.       data airt, airv, airs, airsum/273.11, 1.716e-5, 110.56, 383.67/
  417.       data hfu, wmfu /5.7125e5, 44.09620/
  418.       data hffu, hfco2, hfh2o / -2.35510e6, -8.94114e6, -1.34228e7/
  419.       data gamma/ 1.40 /
  420.       data acpair, bcpair, ccpair/868.3073, 0.2054, -51.846e-6/
  421.       data grav, kgrav/9.80665,0/
  422.       data beta, rdfctr/0.85, 0.1/
  423.       data bta/0.1024/
  424.       data refu, refv, refw, refc/1.0,1.0,1.0,1.0/
  425.       data nitr/1,1,1,1,1/
  426.       data tolr/0.01,0.001,0.01,0.01,0.01/
  427.       data nitke/3/
  428.       data knorth/0/
  429.       data ind1,ind2,ind3,ind4/1,1,1,1/
  430.       data alft, alfrho/0.95,0.95/
  431.       end
  432. C======================================================================C
  433.       subroutine cflux (igrl)
  434. C                                                                      C
  435. C     Purpose:  Calculate values of the mass fluxes and contravariant  C
  436. C               velocities at the cell faces.                          C
  437. C                                                                      C
  438. C======================================================================C
  439.       include 'UIFlow.com'
  440.       include 'UIFlow.indx'
  441. C
  442. C.... Calculate flux perpendicular to eta coordinate.
  443. C
  444.       call grad (igrl, p)
  445.       relxm = 1.0 - relx(1)
  446. C
  447.       do 11 j  = 2, jmax1
  448.          ioff  = (j-1) * imaxl + ibeg
  449.        do 10 i = 2, imax2
  450.            ij  = i + ioff
  451.           ij1  = ij + 1
  452.         cu(ij) = a11(ij) *
  453.      1         ( fx (ij) * ( apu(ij1)+duy(ij1)*dpdy(ij1)/app(ij1) )
  454.      2         + fx1(ij) * ( apu(ij) +duy(ij) *dpdy(ij) /app(ij) ) )
  455.      3         + a12(ij) *
  456.      4         ( fx (ij) * ( apv(ij1)+dvy(ij1)*dpdy(ij1)/ap(ij1))
  457.      5         + fx1(ij) * ( apv(ij) +dvy(ij) *dpdy(ij) /ap(ij)  ) )
  458.      6         + ( a11(ij)*a11(ij)*(fx(ij)/app(ij1) + fx1(ij)/app(ij))
  459.      7         +   a12(ij)*a12(ij)*(fx(ij)/ap(ij1)  + fx1(ij)/ap(ij)))
  460.      8         * ( p(ij) - p(ij1) ) + relxm  * cu(ij)
  461.         cx(ij) = cu(ij) * ( amax1( sign( rho(ij) ,  cu(ij) ), 0.)
  462.      1                  +   amax1( sign( rho(ij1), -cu(ij) ), 0.) )
  463. 10     continue
  464. 11    continue
  465. C
  466. C.... Calculate flux perpendicular to zi coordinate.
  467. C
  468.       do 21 j  = 2, jmax2
  469.          ioff  = (j-1) * imaxl + ibeg
  470.     ioffp  = ioff + imaxl
  471.        do 20 i = 2, imax1
  472.            ij  = i + ioff
  473.       ijp  = i + ioffp
  474.         cv(ij) = a21(ij) *
  475.      1         ( fy (ij) * ( apu(ijp) + dux(ijp)*dpdx(ijp)/app(ijp) )
  476.      2         + fy1(ij) * ( apu(ij)  + dux(ij) *dpdx(ij) /app(ij)) )
  477.      3         + a22(ij) *
  478.      4         ( fy(ij)  * ( apv(ijp) + dvx(ijp)*dpdx(ijp)/ap(ijp) )
  479.      5         + fy1(ij) * ( apv(ij)  + dvx(ij) *dpdx(ij) /ap(ij)  ) )
  480.      6         + ( a21(ij)*a21(ij)*(fy(ij)/app(ijp) + fy1(ij)/app(ij))
  481.      7         +   a22(ij)*a22(ij)*(fy(ij)/ap(ijp)  + fy1(ij)/ap(ij)))
  482.      8         * ( p(ij) - p(ijp) ) + relxm  * cv(ij)
  483.         cy(ij) = cv(ij) * ( amax1( sign( rho(ij) ,  cv(ij) ), 0.)
  484.      1                  +   amax1( sign( rho(ijp), -cv(ij) ), 0.) )
  485. 20     continue
  486. 21    continue
  487.       return
  488.       end
  489. C======================================================================C
  490.       subroutine cfluxc (igrl)
  491. C                                                                      C
  492. C     Purpose:  Make corrections to contravariant velocities and mass  C
  493. C               fluxes.  For use when iterations are being performed   C
  494. C               on the coarse grids.                                   C
  495. C                                                                      C
  496. C======================================================================C
  497.       include 'UIFlow.com'
  498. C
  499.       imaxc  =  imax(igrl)
  500.       jmaxc  =  jmax(igrl)
  501.       imaxf  =  imax(igrl+1)
  502.       jmaxf  =  jmax(igrl+1)
  503.       imaxc1 =  imaxc - 1
  504.       jmaxc1 =  jmaxc - 1
  505.       imaxc2 =  imaxc - 2
  506.       jmaxc2 =  jmaxc - 2
  507.       ibegc  =  nbeg(igrl)
  508. C
  509.       do 6  jc  = 1, jmaxc
  510.        do 5  ic = 1, imaxc
  511.             ijc = ic + (jc-1) * imaxc + ibegc
  512.         dx(ijc) = 0.0
  513.         dy(ijc) = 0.0
  514.         x3(ijc) = 0.0
  515.         su(ijc) = 0.0
  516. 5      continue
  517. 6     continue
  518. C
  519. C.... Evaluate corrections to finer grid.
  520. C
  521.       do 11 jc   = 2, jmaxc1
  522.        do 10 ic  = 2, imaxc1
  523.             ijc  = ic + (jc - 1) * imaxc + ibegc
  524.         ijf  = iru (ijc)
  525.            ijf1  = ijf - 1
  526.            ijfm  = ijf - imaxf
  527.           ijfm1  = ijf - imaxf - 1
  528.         apm(ijc) = ap(ijc)
  529.         dx(ijc)  = u(ijc) - 0.25 * (u(ijf) + u(ijf1) + u(ijfm) +
  530.      1             u(ijfm1))
  531.         dy(ijc)  = v(ijc) - 0.25 * (v(ijf) + v(ijf1) + v(ijfm) +
  532.      1             v(ijfm1))
  533.         x3(ijc)  = p(ijc) - 0.25 * (p(ijf) + p(ijf1) + p(ijfm) +
  534.      1             p(ijfm1))
  535. 10     continue
  536. 11    continue
  537. C
  538.       call bcor (igrl, x3)
  539. C
  540.       call trsrc (igrl, dx)
  541. C
  542.       do 21 j   = 2, jmaxc1
  543.           ioff  = (j-1) * imaxc + ibegc
  544.        do 20 i  = 2, imaxc1
  545.             ij  = i + ioff
  546.         apu(ij) = (ae(ij) * dx(ij+1)      + aw(ij) * dx(ij-1) +
  547.      1             an(ij) * dx(ij+imaxc)  + as(ij) * dx(ij-imaxc) 
  548.      2          +  resux(ij) + su(ij))/app(ij)
  549. 20     continue
  550. 21    continue
  551. C
  552.       call trsrc (igrl, dy)
  553. C
  554.       do 23 j   = 2, jmaxc1
  555.           ioff  = (j-1) * imaxc + ibegc
  556.        do 22 i  = 2, imaxc1
  557.             ij  = i + ioff
  558.         apv(ij) = (ae(ij) * dy(ij+1)      + aw(ij) * dy(ij-1) +
  559.      1             an(ij) * dy(ij+imaxc)  + as(ij) * dy(ij-imaxc) 
  560.      2          +  resvx(ij) + su(ij))/ap(ij)
  561. 22     continue
  562. 23    continue
  563. C
  564. C.... Restrict fluxes again from fine grid.
  565. C
  566.       call restfl (igrl+1)
  567.       call grad (igrl, x3)
  568. C
  569.       do 31 j  = 2, jmaxc1
  570.          ioff  = (j-1) * imaxc + ibegc
  571.         ioffp  = ioff + 1
  572.        do 30 i = 2, imaxc2
  573.            ij  = i + ioff
  574.           ij1  = i + ioffp
  575.         cu(ij) = a11(ij) *
  576.      1         ( fx (ij) * ( apu(ij1)+duy(ij1)*dpdy(ij1)/app(ij1) )
  577.      2         + fx1(ij) * ( apu(ij) +duy(ij) *dpdy(ij) /app(ij) ) )
  578.      3         + a12(ij) *
  579.      4         ( fx (ij) * ( apv(ij1)+dvy(ij1)*dpdy(ij1)/ap(ij1))
  580.      5         + fx1(ij) * ( apv(ij) +dvy(ij) *dpdy(ij) /ap(ij)  ) )
  581.      6         + ( a11(ij)*a11(ij)*(fx(ij)/app(ij1) + fx1(ij)/app(ij))
  582.      7         +   a12(ij)*a12(ij)*(fx(ij)/ap(ij1)  + fx1(ij)/ap(ij)))
  583.      8         * ( x3(ij) - x3(ij1) ) + cu(ij)
  584.         cx(ij) = cu(ij)  * (amax1 (sign (rho(ij) ,  cx(ij)), 0.0)
  585.      1                   +  amax1 (sign (rho(ij1), -cx(ij)), 0.0))
  586. 30     continue
  587. 31    continue
  588. C
  589.       do 41 j  = 2, jmaxc2
  590.           ioff = (j-1) * imaxc + ibegc
  591.      ioffp = ioff + imaxc
  592.        do 40 i = 2, imaxc1
  593.            ij  = i + ioff
  594.       ijp  = i + ioffp
  595.         cv(ij) = a21(ij) *
  596.      1         ( fy (ij) * ( apu(ijp) + dux(ijp)*dpdx(ijp)/app(ijp) )
  597.      2         + fy1(ij) * ( apu(ij)  + dux(ij) *dpdx(ij) /app(ij)) )
  598.      3         + a22(ij) *
  599.      4         ( fy(ij)  * ( apv(ijp) + dvx(ijp)*dpdx(ijp)/ap(ijp) )
  600.      5         + fy1(ij) * ( apv(ij)  + dvx(ij) *dpdx(ij) /ap(ij)  ) )
  601.      6         + ( a21(ij)*a21(ij)*(fy(ij)/app(ijp) + fy1(ij)/app(ij))
  602.      7         +   a22(ij)*a22(ij)*(fy(ij)/ap(ijp)  + fy1(ij)/ap(ij)))
  603.      8         * ( x3(ij) - x3(ijp) ) + cv(ij) 
  604.         cy(ij) = cv(ij)  * (amax1 (sign (rho(ij)  ,  cy(ij)), 0.0)
  605.      1                   +  amax1 (sign (rho(ijp), - cy(ij)), 0.0))
  606. 40     continue
  607. 41    continue
  608.       return
  609.       end    
  610. C======================================================================C
  611.       subroutine cfpmod (igrl)
  612. C                                                                      C
  613. C     Purpose:  Modify coefficients of the pressure correction         C
  614. C               equation to account for boundary condition.            C
  615. C                                                                      C
  616. C======================================================================C
  617.       include 'UIFlow.com'
  618.       include 'UIFlow.indx'
  619. C
  620.       if (kcomp .eq. 1) call pcoefc (igrl)
  621. C
  622.       do 10 j  = 2, jmax1
  623.           ioff = (j-1) * imaxl + ibeg
  624.       ij1  = ioff + 2
  625.       ij2  = ioff + imax1
  626.        ap(ij1) = ap(ij1) - aw(ij1)
  627.        ap(ij2) = ap(ij2) - ae(ij2)
  628.        aw(ij1) = 0.0
  629.        ae(ij2) = 0.0
  630.        dx(ij2) = 0.0
  631. 10    continue
  632. C
  633.       do 20 i  = 2, imax1
  634.           ioff = i + ibeg
  635.       ij1  = ioff + imaxl
  636.       ij2  = ioff + (jmax1-1) * imaxl
  637.        ap(ij1) = ap(ij1) - as(ij1)
  638.        ap(ij2) = ap(ij2) - an(ij2)
  639.        as(ij1) = 0.0
  640.        an(ij2) = 0.0
  641.        dy(ij2) = 0.0
  642. 20    continue
  643. C
  644.       return
  645.       end
  646. C======================================================================C
  647.       subroutine coeff(igrl)
  648. C                                                                      C
  649. C     Purpose:  Calculate coefficients of the discretized equation     C
  650. C               for any flow variable.  The hybrid discretization      C
  651. C               scheme is presently being used.                        C
  652. C                                                                      C
  653. C======================================================================C
  654.       include 'UIFlow.com'
  655.       include 'UIFlow.indx'
  656. C
  657.       do 11 j = 2, jmax1
  658.          ioff = (j-1) * imaxl + ibeg
  659.         ioff1 = ioff - 1
  660.         ioffm = ioff - imaxl
  661.        do 10 i = 2, imax1
  662.            ij  = i + ioff
  663.           ij1  = i + ioff1
  664.           ijm  = i + ioffm
  665.         ae(ij)  = amax1( 0., -cx(ij) , -cx(ij)  * fx (ij)  + dx(ij)  )
  666.         aw(ij)  = amax1( 0.,  cx(ij1),  cx(ij1) * fx1(ij1) + dx(ij1) )
  667.         an(ij)  = amax1( 0., -cy(ij) , -cy(ij)  * fy (ij)  + dy(ij)  )
  668.         as(ij)  = amax1( 0.,  cy(ijm),  cy(ijm) * fy1(ijm) + dy(ijm) )
  669. 10     continue
  670. 11    continue
  671. C
  672.       return
  673.       end
  674. C======================================================================C
  675.       subroutine coeffp (igrl)
  676. C                                                                      C
  677. C     Purpose:  Calculate coefficients of the pressure correction      C
  678. C               equation.                                              C
  679. C                                                                      C
  680. C======================================================================C
  681.       include 'UIFlow.com'
  682.       include 'UIFlow.indx'
  683. C
  684.       do 11 j  = 2, jmax1
  685.           ioff = (j-1) * imaxl + ibeg
  686.          ioff1 = ioff + 1
  687.          ioffp = ioff + imaxl
  688.        do 10 i = 2, imax1
  689.            ij  = i + ioff
  690.           ij1  = i + ioff1
  691.           ijp  = i + ioffp
  692.         ae(ij) = ( a11(ij) * a11(ij)  + a12(ij) * a12(ij) )
  693.      1         * ( fx(ij)  / apm(ij1) + fx1(ij) / apm(ij)  )
  694.         ae(ij) = ae(ij)  * ( amax1( sign( rho(ij)  ,  cx(ij) ), 0.)
  695.      1                   +   amax1( sign( rho(ij1), -cx(ij) ), 0.) )
  696.         dx(ij) = ae(ij)
  697.         an(ij) = ( a21(ij) * a21(ij) + a22(ij) *  a22(ij) )
  698.      1         * ( fy(ij)  / apm(ijp) + fy1(ij) / apm(ij) )
  699.         an(ij) = an(ij)  * ( amax1( sign( rho(ij) ,  cy(ij) ), 0.)
  700.      1                   +   amax1( sign( rho(ijp), -cy(ij) ), 0.) )
  701.         dy(ij) = an(ij)
  702. 10     continue
  703. 11    continue
  704. C
  705.       return
  706.       end
  707. C======================================================================C
  708.       subroutine const
  709. C                                                                      C
  710. C     Purpose:  Compute indices on coarse grids.                       C
  711. C                                                                      C
  712. C======================================================================C
  713.       include 'UIFlow.com'
  714. C
  715.       imax(ngrid) = ncelx + 2
  716.       jmax(ngrid) = ncely + 2
  717.       nbeg (1)    = 0
  718.       if( ngrid .eq. 1 ) return
  719. C
  720.       ngrid1   = ngrid - 1
  721.       do 10 nr = 1, ngrid1
  722.             n  = ngrid - nr
  723.        imax(n) = ncelx / (2 ** nr) + 2
  724.        jmax(n) = ncely / (2 ** nr) + 2
  725. 10    continue
  726. C
  727.       do 15 n  = 2, ngrid
  728.        nbeg(n) = nbeg(n-1) + imax(n-1)*jmax(n-1)
  729. 15    continue
  730. C
  731. C.... Indices for prolongation.
  732. C
  733.       do 25 igr = 1, ngrid1
  734.        imaxc  = imax (igr)
  735.        jmaxc  = jmax (igr)
  736.        ibegc  = nbeg (igr)
  737.        imaxf  = imax(igr+1)
  738.        jmaxf  = jmax(igr+1)
  739.        ibegf  = nbeg(igr+1)
  740.        do 18 j = 1, jmaxc
  741.         do 18 i = 1, imaxc
  742.          ijc     = i + (j-1) * imaxc + ibegc
  743.          ijf     = 2*i-1 + (2*j-2)*imaxf + ibegf
  744.          iru(ijc) = ijf
  745. 18      continue
  746. C
  747. C.... Modify at boundaries.
  748. C
  749.         j = jmaxc
  750.         do 20 i   = 2, imaxc-1
  751.           ijc     = i + (j-1) * imaxc + ibegc
  752.           ijf     = 2*i-1 + (2*j-3)*imaxf + ibegf
  753.          iru(ijc) = ijf
  754. 20      continue
  755.         i = imaxc
  756.         do 21 j   = 2, jmaxc-1
  757.           ijc     = i + (j-1) * imaxc + ibegc
  758.           ijf     = 2*(i-1) + (2*j-2)*imaxf + ibegf
  759.          iru(ijc) = ijf
  760. 21     continue
  761. 25    continue
  762. C
  763. C.... Restrict segment indices.
  764. C
  765.       do 40 nr = 1, ngrid1
  766.             n  = ngrid - nr
  767. C
  768.        do 31 i = 1, nsxm
  769.         jfxm(i,n) = (jfxm (i,n+1)-2)/2 + 2
  770.         jlxm(i,n) = (jlxm (i,n+1)-1)/2 + 1
  771. 31     continue
  772. C
  773.        do 32 i = 1, nsxp
  774.         jfxp(i,n) = (jfxp (i,n+1)-2)/2 + 2
  775.         jlxp(i,n) = (jlxp (i,n+1)-1)/2 + 1
  776. 32     continue
  777. C
  778.        do 33 i = 1, nsym
  779.         ifym(i,n) = (ifym (i,n+1)-2)/2 + 2
  780.         ilym(i,n) = (ilym (i,n+1)-1)/2 + 1
  781. 33     continue
  782. C
  783.        do 34 i = 1, nsyp
  784.         ifyp(i,n) = (ifyp (i,n+1)-2)/2 + 2
  785.         ilyp(i,n) = (ilyp (i,n+1)-1)/2 + 1
  786. 34     continue
  787. C
  788. 40    continue
  789. C
  790.       return
  791.       end
  792. C======================================================================C
  793.       subroutine constr
  794. C                                                                      C
  795. C     Purpose:  Compute the constants required for the combustion      C
  796. C               models.                                                C
  797. C                                                                      C
  798. C======================================================================C
  799.       include 'UIFlow.com'
  800. C
  801. C.... Values for PROPANE ( C3H8 ).
  802. C
  803.       if ( kfuel .eq. 0 ) then
  804.        wmfu   = 44.0962
  805.        wmpr   = 28.3238
  806.        hreact = 46.353e6
  807.        stoic  = 15.5715
  808.        tstoic = 2396.0
  809.        acpfu  = 1420.56
  810.        bcpfu  = 1.7792
  811.        ccpfu  = -3.94184e-4
  812.        acpprd = 882.544
  813.        bcpprd = 1.0537
  814.        ccpprd = -2.22175e-4
  815.        cxx    = 3.0
  816.        hyy    = 8.0
  817.        aa(1)  = 0.10
  818.        bb(1)  = 1.65
  819.        acten(1) = 30000.0
  820.        prexp(1) = 8.6e11
  821.       endif
  822. C
  823. C.... Values for METHANE ( CH4 ).
  824. C
  825.       if ( kfuel .eq. 1 ) then
  826.        wmfu   = 16.04260
  827.        wmpr   = 27.6339
  828.        hreact = 44.164e6
  829.        stoic  = 17.1206
  830.        tstoic = 2329.0
  831.        acpfu  = 1126.59
  832.        bcpfu  = 2.3305
  833.        ccpfu  = -481.2e-6
  834.        acpprd = 891.9962
  835.        bcpprd = 0.3055
  836.        ccpprd = -74.111e-6
  837.        cxx    = 1.0
  838.        hyy    = 4.0
  839.        aa(1)  = -0.3
  840.        bb(1)  = 1.3
  841.        acten(1) = 48400
  842.        prexp(1) = 1.3e8
  843.       endif
  844. C
  845. C.... Values for TOWN GAS ( mixture of methane and ethane ).
  846. C
  847.       if ( kfuel .eq. 2 ) then
  848.        wmfu   = 18.4349
  849.        wmpr   = 27.7509
  850.        hreact = 37.0227e6
  851.        stoic  = 14.376
  852.        tstoic = 2510.0
  853.        acpfu  = 1066.49
  854.        bcpfu  = 1.9408
  855.        ccpfu  = -403.83e-6
  856.        acpprd = 888.335
  857.        bcpprd = 0.3022
  858.        ccpprd = -73.406e-6
  859.        cxx    = 1.0
  860.        hyy    = 4.0
  861.        aa(1)  = -0.3
  862.        bb(1)  = 1.3
  863.        acten(1) = 48400
  864.        prexp(1) = 1.3e8
  865.       endif
  866. C
  867. C.... Calculate constants used for diffusion flame model.
  868. C
  869.       fstoic = 1.0 / ( 1.0 + stoic  )
  870.       rstoic = pref * wmpr / gascon / tstoic
  871.       ysub   = 1.0 / ( 1.0 - fstoic )
  872.       rfuel  = pref * wmfu / gascon / tfuel
  873.       rair   = pref * wmair / gascon / tair
  874.       denom  = fstoic * ( 1.0 - fstoic )
  875.       tprd   = alft * ( tstoic -  ( (1.0 - fstoic) * tair
  876.      1       + fstoic * tfuel ) ) / denom
  877.       rprd   = alfrho * ( rstoic -  ( (1.0 - fstoic) * rair
  878.      1       + fstoic * rfuel ) ) / denom
  879.       phia   = - 1.0 / stoic
  880.       phifma = 1.0 +  ( 1.0 /  stoic )
  881. C
  882.       cpf    = acpfu  + bcpfu * tfuel + ccpfu  * tfuel * tfuel
  883.       cpa    = acpair + bcpair * tair + ccpair * tair  * tair
  884.       enthfu = cpf * tfuel + hreact
  885.       enthox = cpa * tair
  886. C
  887. C.... Calculate constants used for premixed flame model.
  888. C
  889.       rat(1)   = cxx * wmco2 / wmfu
  890.       rat(2)   = ( hyy * wmh2o ) / ( 2.0 * wmfu )
  891.       rat(3)   = cxx * wmox / wmfu + ( hyy * wmox ) / ( 4.0 * wmfu )
  892.       rat(4)   = 3.0 * wmco2 / wmfu
  893.       rat(5)   = 4.0 * wmh2o / wmfu
  894. C
  895.       ab(1)    = aa(1) + bb(1)
  896.       prexp(1) = 1000.0 * wmfu * prexp(1) * ( 0.001 )**ab(1)
  897.      1         * ( 1.0 / wmfu )**aa(1) * ( 1.0 / wmox )**bb(1)
  898.       acten(1) = acten(1) * 4.1868 / 8.3143
  899. C
  900. C....Ensure that initial estimate of mixture fraction is LARGER than
  901. C    the fuel mass fraction.
  902. C
  903.       if ( fugs .gt. fgs ) then
  904.        temp = fgs
  905.        fgs  = fugs
  906.        fugs = temp
  907.       endif
  908. C
  909. C....Specify initial estimates of mass fractions based on the given
  910. C    estimates of mixture fraction and fuel mass fraction.
  911. C
  912.       o2gs  = rox * ( 1.0 - fgs ) - rat(3) * ( fgs - fugs )
  913.       co2gs = rat(1) * ( fgs - fugs )
  914.       h2ogs = rat(2) * ( fgs - fugs )
  915. C
  916.       return
  917.       end
  918. C======================================================================C
  919.       subroutine cpenth (igrl)
  920. C                                                                      C
  921. C     Purpose:  Calculate the specific heat appropriate for            C
  922. C               determination of the fluid temperature.                C
  923. C                                                                      C
  924. C======================================================================C
  925.       include 'UIFlow.com'
  926.       include 'UIFlow.indx'
  927. C
  928. C.... Cp for air only.
  929. C
  930.       if (model .eq. 1) then
  931.        do 11 j   = 1, jmaxl
  932.           ioff   = (j-1) * imaxl + ibeg
  933.         do 10 i  = 1, imaxl
  934.              ij  = i + ioff
  935.          cph(ij) = acpair + bcpair * t(ij) + ccpair * t(ij) * t(ij)
  936. 10      continue
  937. 11     continue
  938. C
  939.       elseif (model .eq. 2) then
  940. C
  941. C.... Cp for mixture of fuel and air.
  942. C
  943.        do 21 j   = 1, jmaxl
  944.            ioff  = (j-1) * imaxl + ibeg
  945.         do 20 i  = 1, imaxl
  946.              ij  = i + ioff
  947.          cpfu    = acpfu  + bcpfu  * t(ij) + ccpfu  * t(ij) * t(ij)
  948.          cpox    = acpox  + bcpox  * t(ij) + ccpox  * t(ij) * t(ij)
  949.          cpn2    = acpn2  + bcpn2  * t(ij) + ccpn2  * t(ij) * t(ij)
  950.          cpco2   = acpco2 + bcpco2 * t(ij) + ccpco2 * t(ij) * t(ij)
  951.          cph2o   = acph2o + bcph2o * t(ij) + ccph2o * t(ij) * t(ij)
  952.          cph(ij) = yfu(ij) * cpfu  +  yo2(ij)  * cpox
  953.      1           + yn2(ij) * cpn2  +  yco2(ij) * cpco2
  954.      2           + yh2o(ij)* cph2o
  955. 20      continue
  956. 21     continue
  957.       endif
  958. C
  959.       return
  960.       end
  961. C======================================================================C
  962.       subroutine cpgrad (igrl)
  963. C                                                                      C
  964. C     Purpose:  Compute the cross derivative terms in the pressure     C
  965. C               correction equation.                                   C
  966. C                                                                      C
  967. C======================================================================C
  968.       include 'UIFlow.com'
  969.       include 'UIFlow.indx'
  970. C
  971.       do 6 j  = 1, jmaxl
  972.         ioff  = (j-1) * imaxl + ibeg
  973.        do 5 i = 1, imaxl
  974.            ij = i + ioff
  975.         x1(ij) = 0.0
  976.         x2(ij) = 0.0
  977. 5      continue
  978. 6     continue
  979. C
  980.       do 11 j  = 2, jmax1
  981.          ioff  = (j-1) * imaxl + ibeg
  982.         ioffp  = ioff + 1
  983.        do 10 i = 2, imax2
  984.            ij  = i + ioff
  985.           ij1  = i + ioffp
  986.         x1(ij) = a11(ij) * (fx(ij) * duy(ij1) * dpdy(ij1)/apm(ij1) +
  987.      1                      fx1(ij)* duy(ij)  * dpdy(ij) /apm(ij)) +
  988.      2           a12(ij) * (fx(ij) * dvy(ij1) * dpdy(ij1)/apm(ij1) +
  989.      3                      fx1(ij)* dvy(ij)  * dpdy(ij) /apm(ij))
  990. 10     continue
  991. 11    continue
  992. C
  993.       do 21 j  = 2, jmax2
  994.          ioff  = (j-1) * imaxl + ibeg
  995.     ioffp  = ioff + imaxl
  996.        do 20 i = 2, imax1
  997.            ij  = i + ioff
  998.       ijp  = i + ioffp
  999.         x2(ij) = a21(ij) * (fy(ij) * dux(ijp) * dpdx(ijp)/apm(ijp) +
  1000.      1                      fy1(ij)* dux(ij)  * dpdx(ij) /apm(ij)) +
  1001.      4           a22(ij) * (fy(ij) * dvx(ijp) * dpdx(ijp)/apm(ijp) +
  1002.      5                      fy1(ij)* dvx(ij)  * dpdx(ij) /apm(ij))
  1003. 20     continue
  1004. 21    continue
  1005. C
  1006.       return
  1007.       end
  1008. C======================================================================C
  1009.       subroutine dens (igrl)
  1010. C                                                                      C
  1011. C     Purpose:  Calculate the fluid density based on the assumption    C
  1012. C               of ideal gas behavior.  Allocation is also made for    C
  1013. C               low mach number flames in which the base pressure of   C
  1014. C               the system determines the thermodynamic state.         C
  1015. C                                                                      C
  1016. C======================================================================C
  1017.       include 'UIFlow.com'
  1018.       include 'UIFlow.indx'
  1019. C
  1020.       relxm = 1.0 - relx(11)
  1021. C
  1022.       if ( kcomp .eq. 0 ) then
  1023. C
  1024. C.... Low mach number case ( Charles' law is applicable ).
  1025. C
  1026.        do 11 j   = 2, jmax1
  1027.           ioff   = (j-1) * imaxl + ibeg
  1028.         do 10 i  = 2, imax1
  1029.              ij  = i + ioff
  1030.          rhl     = pref * wmol(ij) / ( gascon * t(ij) )
  1031.          rho(ij) = relx(11) * rhl + relxm * rho(ij)
  1032. 10      continue
  1033. 11     continue
  1034. C
  1035.       elseif ( kcomp .eq. 1 ) then
  1036. C
  1037. C.... Density determination for a compressible flow.
  1038. C
  1039.        do 21 j   = 2, jmax1
  1040.           ioff   = (j-1) * imaxl + ibeg
  1041.         do 20 i  = 2, imax1
  1042.              ij  = i + ioff
  1043.          rhl     = (p(ij) + pref) * wmol(ij) / ( gascon * t(ij) )
  1044.          rho(ij) = relx(11) * rhl + relxm * rho(ij)
  1045. 20      continue
  1046. 21     continue
  1047.       endif
  1048. C
  1049.       return
  1050.       end
  1051. C======================================================================C
  1052.       subroutine dflux (igrl)
  1053. C                                                                      C
  1054. C     Purpose:  Compute the diffusion contributions to the             C
  1055. C               coefficients in the discretized equations.             C
  1056. C                                                                      C
  1057. C======================================================================C
  1058.       include 'UIFlow.com'
  1059.       include 'UIFlow.indx'
  1060. C
  1061.       do 11 j  = 1, jmax1
  1062.           ioff = (j-1) * imaxl + ibeg
  1063.          ioffp = ioff + imaxl
  1064.        do 10 i = 1, imax1
  1065.            ij  = i + ioff
  1066.           ijp  = i + ioffp
  1067.         dx(ij) = (fx(ij) * gam(ij+1) + fx1(ij) * gam(ij)) * q11(ij)
  1068.         dy(ij) = (fy(ij) * gam(ijp)  + fy1(ij) * gam(ij)) * q22(ij)
  1069. 10     continue
  1070. 11    continue
  1071. C
  1072.       return
  1073.       end
  1074. C======================================================================C
  1075.       subroutine enth (igrl)
  1076. C                                                                      C
  1077. C     Purpose:  Calculate initial sensible enthalpy field.             C
  1078. C                                                                      C
  1079. C======================================================================C
  1080.       include 'UIFlow.com'
  1081.       include 'UIFlow.indx'
  1082. C
  1083.       if (model .eq. 1) then
  1084. C
  1085. C.... For air only.
  1086. C
  1087.       hmix = rox * hox  +  (1.0 - rox) * hn2
  1088. C
  1089.        do 11 j  = 1, jmaxl
  1090.           ioff  = (j-1) * imaxl + ibeg
  1091.         do 10 i = 1, imaxl
  1092.              ij = i + ioff
  1093.          h(ij)  = cph(ij) * t(ij) - hmix
  1094. 10      continue
  1095. 11     continue
  1096. C
  1097.       elseif (model .eq. 2) then
  1098. C
  1099. C.... For constituents of the one step chemical reaction.
  1100. C
  1101.        do 21 j  = 1, jmaxl
  1102.           ioff  = (j-1) * imaxl + ibeg
  1103.         do 20 i = 1, imaxl
  1104.              ij = i + ioff
  1105.          hmix   =  yfu(ij)  * hfu  +  yo2(ij)  * hox 
  1106.      1          +  yn2(ij)  * hn2  +  yco2(ij) * hco2
  1107.      2          +  yh2o(ij) * hh2o
  1108.          h(ij)  =  cph(ij)  * t(ij) - hmix
  1109. 20      continue
  1110. 21     continue
  1111.       endif
  1112. C
  1113.       return
  1114.       end
  1115. C======================================================================C
  1116.       subroutine extrpl(igrl)
  1117. C                                                                      C
  1118. C     Purpose:  Enforce boundary conditions by extrapolating interior  C
  1119. C               values.                                                C
  1120. C                                                                      C
  1121. C======================================================================C
  1122.       include 'UIFlow.com'
  1123. C
  1124.       call xmextr (igrl)
  1125.       call xpextr (igrl)
  1126.       call ymextr (igrl)
  1127.       call ypextr (igrl)
  1128.       if ( kadj .eq. 0 ) call mssout (igrl)
  1129. C
  1130.       return
  1131.       end
  1132. C======================================================================C
  1133.       subroutine fstchm (igrl)
  1134. C                                                                      C
  1135. C     Purpose:  Calculate quantities pertinent to the diffusion flame  C
  1136. C               model.                                                 C
  1137. C                                                                      C
  1138. C======================================================================C
  1139.       include 'UIFlow.com'
  1140.       include 'UIFlow.indx'
  1141. C
  1142. C.... Properties for a LAMINAR reacting flow assuming a mixing limited
  1143. C     reaction.
  1144. C
  1145.       call scalar (igrl, 8, f)
  1146. C
  1147.       relxm  = 1.0 - relx(5)
  1148.       rwmfu  = 1.0 / wmfu
  1149.       rwmair = 1.0 / wmair
  1150.       rwmpr  = 1.0 / wmpr
  1151. C
  1152.       do 11 j     = 2, jmax1
  1153.           ioff    = (j-1) * imaxl + ibeg
  1154.        do 10 i    = 2, imax1
  1155.              ij   = i + ioff
  1156.         phi       = f(ij) * phifma + phia
  1157.         yfu  (ij) = amax1( 0.0, phi )
  1158.         yo2  (ij) = ( yfu(ij) - phi ) * stoic
  1159.         yco2 (ij) = 1.0 - yfu(ij) - yo2(ij)
  1160.         h    (ij) = f(ij) * enthfu + (1.0 - f(ij)) * enthox
  1161.         cpf       = acpfu  + bcpfu  * t(ij) + ccpfu  * t(ij) * t(ij)
  1162.         cpa       = acpair + bcpair * t(ij) + ccpair * t(ij) * t(ij)
  1163.         cprd      = acpprd + bcpprd * t(ij) + ccpprd * t(ij) * t(ij)
  1164.         cph  (ij) = yfu(ij)  * cpf + yo2(ij) * cpa
  1165.      1            + yco2(ij) * cprd
  1166.         tl        = ( h(ij) - yfu(ij) * hreact ) / cph(ij)
  1167.         t(ij)     = relx(5) * tl + relxm * t(ij)
  1168.         rwmol     = yfu(ij)*rwmfu + yo2(ij)*rwmair + yco2(ij)*rwmpr
  1169.         wmol(ij)  = 1.0 / rwmol
  1170. 10     continue
  1171. 11    continue
  1172.       call dens (igrl)
  1173. C
  1174.       if (klam .eq. 0) then
  1175. C
  1176. C.... Properties for a TURBULENT reacting flow assuming a mixing
  1177. C     limited reaction.
  1178. C
  1179.        call scalar (igrl, 10, g)
  1180. C
  1181.        relxmr = 1.0 - relx(11)
  1182.        relxmt = 1.0 - relx(5)
  1183.        stcp1  = 1.0 + stoic
  1184. C
  1185.        do 31 j  = 2, jmax1
  1186.           ioff  = (j-1) * imaxl + ibeg
  1187.         do 30 i = 2, imax1
  1188.             ij  = i + ioff
  1189.          fprme    =  sqrt( g(ij) )
  1190.          zsubs    = abs( fstoic - f(ij) ) / fprme
  1191.          unmxd    = 0.45 * exp ( - zsubs )
  1192.          unprd    = ysub * fprme * unmxd
  1193.          yco2(ij) = yco2(ij) - stcp1 * unprd
  1194.          if ( yco2(ij) .lt. 0 ) then
  1195.           unmxd = 0.5 * unmxd
  1196.           unprd = 0.5 * unprd
  1197.          endif
  1198.          yfu(ij)  = yfu(ij) + unprd
  1199.          yo2(ij)  = yo2(ij) + stoic * unprd
  1200.          yco2(ij) = 1.0 - yfu(ij) - yo2(ij)
  1201.          tl       = t(ij) - tprd * fprme * unmxd
  1202.          t(ij)    = relx(5) * tl + relxmt * t(ij)
  1203.          rrho     = 1.0 / rho(ij) - rprd * fprme * unmxd
  1204.          rhl      = 1.0 / rrho
  1205.          rho(ij)  = relx(11) * rhl + relxmr * rho(ij)
  1206.          rwmol    = yfu(ij)*rwmfu + yo2(ij)*rwmair + yco2(ij)*rwmpr
  1207.          wmol(ij) = 1.0 / rwmol
  1208. 30      continue
  1209. 31     continue
  1210.       endif
  1211. C
  1212.       return
  1213.       end
  1214. C======================================================================C
  1215.       subroutine ftout(i      jfl = jfxp(n,igrl)
  1216.       jll = jlxp(n,igrl)
  1217.        call  wallg (igrl,i,i,jfl,jll,1,xxic)
  1218. 20    continue
  1219. C
  1220.       j = 2
  1221.       do 30 n = 1, nsym
  1222.        if ( kbym(n) .ne. 1 ) goto 30
  1223.          ifl  = ifym(n,igrl)
  1224.      ill  = ilym(n,igrl)
  1225.        call  wallg (igrl,ifl,ill,j,j,-imaxl,yetac)
  1226. 30    continue
  1227. C
  1228.       j = jmax1
  1229.       do 40 n = 1, nsyp
  1230.        if ( kbyp(n) .ne. 1 ) goto 40
  1231.          ifl  = ifyp(n,igrl)
  1232.      ill  = ilyp(n,igrl)
  1233.        call  wallg (igrl,ifl,ill,j,j,imaxl,yetac)
  1234. 40    continue
  1235. C
  1236.       return
  1237.       end
  1238. C======================================================================C
  1239.       subroutine geomc1 (igrl)
  1240. C                                                                      C
  1241. C     Purpose:  Calculate geometry terms associated with the           C
  1242. C               transformation to general curvilinear coordinates.     C
  1243. C               Geometry terms stored at cell centers are              C
  1244. C               calculated here.                                       C
  1245. C                                                                      C
  1246. C======================================================================C
  1247.       include 'UIFlow.com'
  1248.       include 'UIFlow.indx'
  1249. C
  1250.       do 11 j   = 2, jmax1
  1251.          ioff   = (j-1) * imaxl + ibeg
  1252.         ioff1   = ioff - imaxl
  1253.        do 10 i  = 2, imax1
  1254.             ij  = i + ioff
  1255.            ij1  = i + ioff1
  1256.         xzi     = 0.5 * ( x(ij)  + x(ij1)  -  x(ij-1) - x(ij1-1) )
  1257.         yzi     = 0.5 * ( y(ij)  + y(ij1)  -  y(ij-1) - y(ij1-1) )
  1258.         xeta    = 0.5 * ( x(ij)  + x(ij-1) -  x(ij1)  - x(ij1-1) )
  1259.         yeta    = 0.5 * ( y(ij)  + y(ij-1) -  y(ij1)  - y(ij1-1) )
  1260.         ajb(ij) =   xzi * yeta - xeta * yzi
  1261.         q12(ij) = - yzi * yeta - xzi  * xeta
  1262.         xxic (ij) = sqrt( xzi  * xzi  + yzi  * yzi  )
  1263.         yetac(ij) = sqrt( yeta * yeta + xeta * xeta )
  1264.         dux(ij) =   yeta
  1265.         duy(ij) = - yzi
  1266.         dvx(ij) = - xeta
  1267.         dvy(ij) =   xzi
  1268.         q12(ij) =   q12(ij) / ajb(ij)
  1269.         q21(ij) =   q12(ij)
  1270. 10     continue
  1271. 11    continue
  1272. C
  1273. C.... Change geometry terms if the grid describes an AXISYMMETRIC domain.
  1274. C
  1275.       if ( kplax .eq. 1 ) return
  1276. C
  1277.       do 21 j   = 2, jmax1
  1278.          ioff   = (j-1) * imaxl + ibeg
  1279.         ioff1   = ioff - imaxl
  1280.        do 20 i  = 2, imax1
  1281.             ij  = i + ioff
  1282.            ij1  = i + ioff1
  1283.         rav     = 0.25 * ( r(ij) + r(ij-1) + r(ij1) + r(ij1-1) )
  1284.         dux(ij) = dux(ij) * rav
  1285.         duy(ij) = duy(ij) * rav
  1286.         dvx(ij) = dvx(ij) * rav
  1287.         dvy(ij) = dvy(ij) * rav
  1288.         q12(ij) = q12(ij) * rav
  1289.         q21(ij) = q12(ij)
  1290. 20     continue
  1291. 21    continue
  1292. C
  1293.       return
  1294.       end
  1295. C======================================================================C
  1296.       subroutine geomc2 (igrl)
  1297. C                                                                      C
  1298. C     Purpose:  Modify the Jacobian for axisymmetric systems.          C
  1299. C               This modification must NOT be done in Geomc1.          C
  1300. C                                                                      C
  1301. C======================================================================C
  1302.       include 'UIFlow.com'
  1303.       include 'UIFlow.indx'
  1304. C
  1305.       if (kplax .eq. 1) return
  1306. C
  1307.       do 11 j   = 2, jmax1
  1308.           ioff  = (j-1) * imaxl + ibeg
  1309.          ioff1  = ioff - imaxl
  1310.        do 10 i  = 2, imax1
  1311.           ij    = i + ioff
  1312.           ij1   = i + ioff1
  1313.         rav     =  0.25 * ( r(ij) + r(ij-1) + r(ij1) + r(ij1-1) )
  1314.         ajb(ij) =  ajb(ij) * rav
  1315. 10     continue
  1316. 11    continue
  1317. C
  1318.       return
  1319.       end
  1320. C======================================================================C
  1321.       subroutine geomm
  1322. C                                                                      C
  1323. C     Purpose:  Compute all transformation metrics and interpolation   C
  1324. C               factors.                                               C
  1325. C                                                                      C
  1326. C======================================================================C
  1327.       include 'UIFlow.com'
  1328. C
  1329.       do 200 igr = 1,ngrid 
  1330.        call geomc1 (igr)
  1331.        call geomx  (igr)
  1332.        call geomy  (igr)
  1333.        call geomc2 (igr)
  1334.        call fxx    (igr)
  1335.        call fyy    (igr)
  1336. 200   continue
  1337. C
  1338.       return
  1339.       end
  1340. C======================================================================C
  1341.       subroutine geomx (igrl)
  1342. C                                                                      C
  1343. C     Purpose:  Calculate geometry terms associated with the           C
  1344. C               transformation to general curvilinear coordinates.     C
  1345. C               Geometry terms stored at cell faces of constant zi     C
  1346. C               are calculated here.                                   C
  1347. C                                                                      C
  1348. C======================================================================C
  1349.       include 'UIFlow.com'
  1350.       include 'UIFlow.indx'
  1351. C
  1352.       do 11 j   = 2, jmax1
  1353.          ioff   = (j-1) * imaxl + ibeg
  1354.         ioff1   = ioff - imaxl
  1355.        do 10 i  = 1, imax1
  1356.             ij  = i + ioff
  1357.            ij1  = i + ioff1
  1358.         xeta    = ( x(ij) - x(ij1) )
  1359.         yeta    = ( y(ij) - y(ij1) )
  1360.         a11(ij) =   yeta
  1361.         a12(ij) = - xeta
  1362.         q11(ij) = a11(ij) * a11(ij)  +  a12(ij) * a12(ij)
  1363.         q11(ij) = 2.0 * q11(ij) / ( ajb(ij) + ajb(ij+1) )
  1364. 10     continue
  1365. 11    continue
  1366. C
  1367. C.... Change geometry terms if the grid describes an AXISYMMETRIC domain.
  1368. C
  1369.       if ( kplax .eq. 1 ) return
  1370. C
  1371.       do 21 j   = 2, jmax1
  1372.          ioff   = (j-1) * imaxl + ibeg
  1373.         ioff1   = ioff - imaxl
  1374.        do 20 i  = 1, imax1
  1375.             ij  = i + ioff
  1376.            ij1  = i + ioff1
  1377.         rav     = 0.5 * ( r(ij) + r(ij1) )
  1378.         a11(ij) = a11(ij) * rav
  1379.         a12(ij) = a12(ij) * rav
  1380.         q11(ij) = q11(ij) * rav
  1381. 20     continue
  1382. 21    continue
  1383. C
  1384.       return
  1385.       end
  1386. C======================================================================C
  1387.       subroutine geomy (igrl)
  1388. C                                                                      C
  1389. C     Purpose:  Calculate geometry terms associated with the           C
  1390. C               transformation to general curvilinear coordinates.     C
  1391. C               Geometry terms stored at cell faces of constant eta    C
  1392. C               are calculated here.                                   C
  1393. C                                                                      C
  1394. C======================================================================C
  1395.       include 'UIFlow.com'
  1396.       include 'UIFlow.indx'
  1397. C
  1398.       do 11 j   = 1, jmax1
  1399.          ioff   = (j-1) * imaxl + ibeg
  1400.        do 10 i  = 2, imax1
  1401.             ij  = i + ioff
  1402.         xzi     = ( x(ij) - x(ij-1) )
  1403.         yzi     = ( y(ij) - y(ij-1) )
  1404.         a21(ij) = - yzi
  1405.         a22(ij) =   xzi
  1406.         q22(ij) = a21(ij) * a21(ij)  +  a22(ij) * a22(ij)
  1407.         q22(ij) = 2.0 * q22(ij) / ( ajb(ij) + ajb(ij+imaxl) )
  1408. 10     continue
  1409. 11    continue
  1410. C
  1411. C.... Change geometry terms if the grid describes an AXISYMMETRIC domain.
  1412. C
  1413.       if ( kplax .eq. 1 ) return
  1414. C
  1415.       do 21 j   = 1, jmax1
  1416.          ioff   = (j-1) * imaxl + ibeg
  1417.        do 20 i  = 2, imax1
  1418.             ij  = i + ioff
  1419.         rav     = 0.5 * ( r(ij) + r(ij-1) )
  1420.         a21(ij) = a21(ij) * rav
  1421.         a22(ij) = a22(ij) * rav
  1422.         q22(ij) = q22(ij) * rav
  1423. 20     continue
  1424. 21    continue
  1425.       return
  1426.       end
  1427. C======================================================================C
  1428.       subroutine grad (igrl,q)
  1429. C                                                                      C
  1430. C     Purpose:  Calculate the NEGATIVE of the zi and eta gradients     C
  1431. C               for the given variable.                                C
  1432. C                                                                      C
  1433. C======================================================================C
  1434.       include 'UIFlow.com'
  1435.       dimension q(*)
  1436.       include 'UIFlow.indx'
  1437. C
  1438.       do 12 j   = 2, jmax1
  1439.           ioff  = (j-1) * imaxl + ibeg
  1440.        do 11 i   = 2, imax1
  1441.             ij   = i + ioff
  1442.            ij1   = ij - 1
  1443.            ijm   = ij - imaxl
  1444.         dpdx(ij) = fx(ij1) * q(ij)   + fx1(ij1) * q(ij1)
  1445.      1           - fx(ij)  * q(ij+1) - fx1(ij)  * q(ij)
  1446.         dpdy(ij) = fy(ijm) * q(ij)   + fy1(ijm) * q(ijm)
  1447.      1           - fy1(ij) * q(ij)   - fy (ij)  * q(ij+imaxl)
  1448. 11     continue
  1449. 12    continue
  1450. C
  1451.       return
  1452.       end
  1453. C======================================================================C
  1454.       subroutine grid
  1455. C                                                                      C
  1456. C     Purpose:  Read in the x,y locations of the finest grid and then  C
  1457. C               calculate the coarser grids.                           C
  1458. C                                                                      C
  1459. C======================================================================C
  1460.       include 'UIFlow.com'
  1461. C
  1462.       imaxl  =  imax  (ngrid)
  1463.       jmaxl  =  jmax  (ngrid)
  1464.       ibeg   =  nbeg  (ngrid)
  1465.       imax1  =  imaxl - 1
  1466.       jmax1  =  jmaxl - 1
  1467. C
  1468.       do 10 i  = 1, imax1
  1469.        do 11 j = 1, jmax1
  1470.             ij = i + (j-1) * imaxl + ibeg
  1471.         read (11,*) x(ij), y(ij)
  1472. 11     continue
  1473. 10    continue
  1474. C
  1475.       do 20 i= 1, imax1
  1476.        do 21 j = 1, jmax1
  1477.            ij  = i + (j-1) * imaxl + ibeg
  1478.         r(ij)  = y(ij) + rin
  1479. 21     continue
  1480. 20    continue
  1481. C
  1482. C.... Calculate coarse grids based on the user prescribed finest grid.
  1483. C
  1484.       ngrid1 = ngrid - 1
  1485.       if (ngrid .eq. 1) return
  1486.       do 50 n = 1, ngrid1
  1487.          nrev    = ngrid - n
  1488.          imaxc   = imax (nrev)
  1489.          jmaxc   = jmax (nrev)
  1490.          imaxf   = imax (nrev + 1)
  1491.          jmaxf   = jmax (nrev + 1)
  1492.          imaxc1  = imaxc - 1
  1493.          jmaxc1  = jmaxc - 1
  1494.          ibegc   = nbeg (nrev)
  1495.          ibegf   = nbeg (nrev + 1)
  1496.        do 30 ic  = 1, imaxc1
  1497.         do 31 jc = 1, jmaxc1
  1498.              if  = 2 * ic - 1
  1499.          jf  = 2 * jc - 1
  1500.             ijf  = if + (jf - 1) * imaxf + ibegf 
  1501.         ijc  = ic + (jc - 1) * imaxc + ibegc 
  1502.          x(ijc)  = x(ijf)
  1503.          y(ijc)  = y(ijf)
  1504.          r(ijc)  = r(ijf)
  1505. 31      continue
  1506. 30     continue
  1507. 50    continue
  1508. C
  1509.       return
  1510.       end
  1511. C======================================================================C
  1512.       subroutine gridp
  1513. C                                                                      C
  1514. C     Purpose:  Read in x,y locations of the finest grid and calculate C
  1515. C               the x,y locations for the coarser grids.  This routine C
  1516. C               is used in conjunction with a grid generated by the    C
  1517. C               pre-processor.                                         C
  1518. C                                                                      C
  1519. C======================================================================C
  1520.       include 'UIFlow.com'
  1521. C
  1522.       imaxl  =  imax  (1)
  1523.       jmaxl  =  jmax  (1)
  1524.       ibeg   =  nbeg  (1)
  1525.       imax1  =  imaxl - 1
  1526.       jmax1  =  jmaxl - 1
  1527. C
  1528.       do 10 i  = 1, imax1
  1529.        do 11 j = 1, jmax1
  1530.            ij  = i + (j-1) * imaxl + ibeg
  1531.         read (11,*) x(ij), y(ij)
  1532. 11     continue
  1533. 10    continue
  1534. C
  1535.       do 20 i  = 1, imax1
  1536.        do 21 j = 1, jmax1
  1537.            ij  = i + (j-1) * imaxl + ibeg
  1538.         r(ij)  = y(ij) + rin
  1539. 21     continue
  1540. 20    continue
  1541. C
  1542. C.... Calculate coarse grids based on the user prescribed finest grid.
  1543. C
  1544.       if (ngrid .eq. 1) return
  1545.       do 50 n = 2, ngrid
  1546.        imaxc  = imax (n-1)
  1547.        jmaxc  = jmax (n-1)
  1548.        imaxf  = imax (n)
  1549.        jmaxf  = jmax (n)
  1550.        imaxc1 = imaxc - 1
  1551.        jmaxc1 = jmaxc - 1
  1552.        ibegc  = nbeg (n-1)
  1553.        ibegf  = nbeg (n)
  1554.        do 30 ic  = 1, imaxc1
  1555.         do 31 jc = 1, jmaxc1
  1556.              if  = 2 * ic - 1
  1557.          jf  = 2 * jc - 1
  1558.             ijf  = if + (jf - 1) * imaxf + ibegf
  1559.         ijc  = ic + (jc - 1) * imaxc + ibegc
  1560. C
  1561.          x(ijf)   = x(ijc)
  1562.          y(ijf)   = y(ijc)
  1563.          r(ijf)   = r(ijc)
  1564. C
  1565.          x(ijf+1) = 0.5*( x(ijc) + x(ijc+1) )
  1566.          y(ijf+1) = 0.5*( y(ijc) + y(ijc+1) )
  1567.          r(ijf+1) = 0.5*( r(ijc) + r(ijc+1) )
  1568. C
  1569.          x(ijf+imaxf) = 0.5*( x(ijc) + x(ijc+imaxc) )
  1570.          y(ijf+imaxf) = 0.5*( y(ijc) + y(ijc+imaxc) )
  1571.          r(ijf+imaxf) = 0.5*( r(ijc) + r(ijc+imaxc) )
  1572. C
  1573.          x(ijf+imaxf+1) = 0.25*( x(ijc)+x(ijc+1)+x(ijc+imaxc)+
  1574.      1                           x(ijc+1+imaxc) )
  1575.          y(ijf+imaxf+1) = 0.25*( y(ijc)+y(ijc+1)+y(ijc+imaxc)+
  1576.      1                           y(ijc+1+imaxc) )
  1577.          r(ijf+imaxf+1) = 0.25*( r(ijc)+r(ijc+1)+r(ijc+imaxc)+
  1578.      1                           r(ijc+1+imaxc) )
  1579. 31      continue
  1580. 30     continue
  1581. 50    continue
  1582. C
  1583.       return
  1584.       end
  1585. C======================================================================C
  1586.       subroutine header
  1587. C                                                                      C
  1588. C     Purpose:  Print a title for program output.                      C
  1589. C                                                                      C
  1590. C======================================================================C
  1591.        write(8,*) 
  1592.        write(8,*) ('*',i = 1,86)
  1593.        write(8,*) ('*',i = 1,86)
  1594.        do 10 j = 1,22
  1595.          write(8,*) '**',(' ',i = 1,82),'**'
  1596.  10    continue
  1597.        write(8,*) '**',(' ',i=1,34),'UIFLOW - 2D',(' ',i=1,34),'**'
  1598.        write(8,*) '**',(' ',i= 1,82),'**' 
  1599.        write(8,*) '**',(' ',i=1,35),'VERSION 2.0 ',(' ',i=1,35),'**' 
  1600.        write(8,*) '**',(' ',i= 1,82),'**' 
  1601.        write(8,*) '**',(' ',i= 1,82),'**' 
  1602.        write(8,*) '**',(' ',i= 1,82),'**' 
  1603.        write(8,*) '**',(' ',i=1,19),'UNIVERSITY of ILLINOIS at
  1604.      1 URBANA - CHAMPAIGN',(' ',i=1,19),'**'
  1605.        write(8,*) '**',(' ',i= 1,82),'**' 
  1606.        write(8,*) '**',(' ',i=1,15),'DEPARTMENT of MECHANICAL and 
  1607.      1 INDUSTRIAL ENGINEERING',(' ',i=1,15),'**'
  1608.        write(8,*) '**',(' ',i= 1,82),'**'
  1609.        write(8,*) '**',(' ',i= 1,82),'**'
  1610.        write(8,*) '**',(' ',i=1,29),'CONTACT - Dr. S.P. Vanka',
  1611.      1 (' ',i=1,29),'**'
  1612.        write(8,*) '**',(' ',i=1,33),'(217) 244 - 8388',
  1613.      1 (' ',i=1,33),'**' 
  1614.        do 20 j = 1,23
  1615.          write(8,*) '**',(' ',i= 1,82),'**'
  1616.  20    continue
  1617.        write(8,*) ('*',i = 1,86)
  1618.        write(8,*) ('*',i = 1,86)
  1619.        write(8,*)
  1620.        return
  1621.        end
  1622. C======================================================================C
  1623.       subroutine init
  1624. C                                                                      C
  1625. C     Purpose:  Specify initial estimate of flowfield variables on     C
  1626. C               all grids.  Also specify inlet conditions on all       C
  1627. C               grids.                                                 C
  1628. C                                                                      C
  1629. C======================================================================C
  1630.       include 'UIFlow.com'
  1631. C
  1632. C.... Specify inlet conditions and initial estimates on the finest grid.
  1633. C
  1634.       call  intern
  1635.       call  lbndry
  1636.       call  rbndry
  1637.       call  tbndry
  1638.       call  bbndry
  1639.       call  mssin
  1640. C
  1641.       call  cpenth (ngrid)
  1642.       call  enth   (ngrid)
  1643. C
  1644. C.... Restrict flow variables to coarser grids.
  1645. C
  1646.       if( ngrid .gt. 1 ) then
  1647.        ngrid1   = ngrid - 1
  1648.        do 10 nr  = 1, ngrid1
  1649.              igr = ngrid - nr + 1 
  1650.         call restv (igr)
  1651.         call restfl(igr)
  1652. C
  1653.         if (klam .eq. 0) then
  1654.          call rests (igr, tke)
  1655.          call rests (igr, tde)
  1656.         endif
  1657. C
  1658.         if (model.gt. 0) then
  1659.          call rests (igr, h)
  1660.          call rests (igr, t)
  1661.          call rests (igr, wmol)
  1662.         endif
  1663. C
  1664.         if (model .eq. 2) then
  1665.          call rests (igr, f)
  1666.          call rests (igr, yfu)
  1667.          call rests (igr, yo2)
  1668.          call rests (igr, yco2)
  1669.          call rests (igr, yh2o)
  1670.          call rests (igr, yn2)
  1671.         endif
  1672. C
  1673.         if (model .eq. 3) then
  1674.          call rests (igr, f)
  1675.          if (klam .eq. 0) call rests (igr, g)
  1676.          call rests (igr, yfu)
  1677.          call rests (igr, yo2)
  1678.          call rests (igr, yco2)
  1679.         endif
  1680. C
  1681.         call extrpl(igr-1)
  1682. 10     continue
  1683.       endif
  1684.       return
  1685.       end
  1686. C======================================================================C
  1687.       subroutine input
  1688. C                                                                      C
  1689. C     Purpose:  Read in all data pertinent to describing the flow      C
  1690. C               to be solved.                                          C
  1691. C                                                                      C
  1692. C======================================================================C
  1693.       include 'UIFlow.com'
  1694. C
  1695. C.... Read general problem data.
  1696. C
  1697.       read (5, *) klam , kcomp , kswrl, kpgrid, model
  1698.       read (5, *) kfuel, knorth, kplax, kadj
  1699.       read (5, *) ngrid, ncelx, ncely
  1700. C
  1701. C.... Read segment data.
  1702. C
  1703.       read (5, *) nsxm
  1704.       do 10 n = 1,nsxm
  1705.        read (5, *) kbxm(n), jfc, jlc
  1706.        jfxm (n, ngrid) = jfc + 1 
  1707.        jlxm (n, ngrid) = jlc + 1 
  1708.        read (5, *) ubxm(n), vbxm(n) , wbxm(n) , dvar   , txm(n)  
  1709.        read (5, *) rhxm(n), fxm(n)  , gxm(n)  , tkxm(n), tdxm(n) 
  1710.        read (5, *) fuxm(n), co2xm(n), h2oxm(n), o2xm(n), wmxm(n)
  1711. 10    continue
  1712. C
  1713.       read (5, *) nsxp
  1714.       do 20 n = 1,nsxp
  1715.        read (5, *) kbxp (n), jfc, jlc
  1716.        jfxp (n, ngrid) = jfc + 1 
  1717.        jlxp (n, ngrid) = jlc + 1 
  1718.        read (5, *) ubxp(n), vbxp(n) , wbxp(n) , dvar   , txp(n)       
  1719.        read (5, *) rhxp(n), fxp(n)  , gxp(n)  , tkxp(n), tdxp(n)
  1720.        read (5, *) fuxp(n), co2xp(n), h2oxp(n), o2xp(n), wmxp(n)
  1721. 20    continue
  1722. C
  1723.       read (5, *) nsym     
  1724.       do 30 n = 1,nsym
  1725.        read (5, *) kbym (n), ifc, ilc
  1726.        ifym (n, ngrid) = ifc + 1 
  1727.        ilym (n, ngrid) = ilc + 1 
  1728.        read (5, *) ubym(n), vbym(n) , wbym(n) , dvar   , tym(n)
  1729.        read (5, *) rhym(n), fym(n)  , gym(n)  , tkym(n), tdym(n)
  1730.        read (5, *) fuym(n), co2ym(n), h2oym(n), o2ym(n), wmym(n)
  1731. 30    continue
  1732. C
  1733.       read (5, *) nsyp
  1734.       do 40 n = 1,nsyp
  1735.        read (5, *) kbyp(n), ifc, ilc
  1736.        ifyp (n, ngrid) = ifc + 1 
  1737.        ilyp (n, ngrid) = ilc + 1 
  1738.        read (5, *) ubyp(n), vbyp(n) , wbyp(n) , dvar   , typ(n)
  1739.        read (5, *) rhyp(n), fyp(n)  , gyp(n)  , tkyp(n), tdyp(n)
  1740.        read (5, *) fuyp(n), co2yp(n), h2oyp(n), o2yp(n), wmyp(n)
  1741. 40    continue
  1742. C
  1743.       read (5, *) ugs, vgs, wgs, rhgs
  1744.       read (5, *) tgs, tkgs, tdgs, fgs, ggs, fugs
  1745.       read (5, *) tfuel, tair
  1746.       read (5, *) (prl(nv)  , nv=1, 11)
  1747.       read (5, *) (prt(nv)  , nv=1, 11)
  1748.       read (5, *) (relx(nv) , nv=1, 11)
  1749.       read (5, *) (nswp(nv) , nv=1, 11)
  1750.       read (5, *) (iprint(i), i=1 , 12)
  1751.       read (5, *) pref, vscty
  1752.       read (5, *) maxitn, tolr(ngrid)
  1753.       read (5, *) rin
  1754.       read (5, *) nxbaf, nybaf, nobs
  1755.       relx(2) =  relx(1)
  1756. C
  1757. C      call xwcommon
  1758. 1001  format (80a1)
  1759. C
  1760.       return
  1761.       end
  1762. C======================================================================C
  1763.       subroutine intern
  1764. C                                                                      C
  1765. C     Purpose:  Initialize internal values of flow variables.          C
  1766. C                                                                      C
  1767. C======================================================================C
  1768.       include 'UIFlow.com'
  1769.       igrl = ngrid
  1770.       include 'UIFlow.indx'
  1771. C
  1772.       if(model .gt. 0) rhgs = pref * wmair / ( gascon * tgs )
  1773. C
  1774.       do 11 i     =  2, imax1
  1775.        do 10 j    =  2, jmax1
  1776.            ijf    =  i + (j-1) * imaxl + ibeg
  1777.         u(ijf)    =  ugs
  1778.         v(ijf)    =  vgs
  1779.         w(ijf)    =  wgs
  1780.         p(ijf)    =  0.0
  1781.         t(ijf)    =  tgs
  1782.         gam(ijf)  =  vscty
  1783.         amu(ijf)  =  vscty
  1784.         amut(ijf) =  cd * rhgs * tkgs * tkgs / ( tdgs + 1.0e-30 )
  1785.         rho(ijf)  =  rhgs
  1786.         tke(ijf)  =  tkgs
  1787.         tde(ijf)  =  tdgs
  1788.         f  (ijf)  =  fgs
  1789.         g  (ijf)  =  ggs
  1790.         yfu(ijf)  =  fugs
  1791.         yco2(ijf) =  co2gs
  1792.         yh2o(ijf) =  h2ogs
  1793.         yo2(ijf)  =  o2gs
  1794.         yn2(ijf)  =  1.0 - fugs - co2gs - h2ogs - o2gs
  1795.         wmol(ijf) =  wmair
  1796.         cu(ijf)   =  a11(ijf) * u(ijf) + a12(ijf) * v(ijf)
  1797.         cv(ijf)   =  a21(ijf) * u(ijf) + a22(ijf) * v(ijf)
  1798.         cx(ijf)   =  rho(ijf) * cu(ijf)
  1799.         cy(ijf)   =  rho(ijf) * cv(ijf)
  1800. 10     continue
  1801. 11    continue
  1802. C
  1803.       return
  1804.       end
  1805. C======================================================================C
  1806.       subroutine lamvis (igrl)
  1807. C                                                                      C
  1808. C     Purpose:  Compute the laminar viscosity based on the             C
  1809. C               Sutherland's law constants for AIR only.               C
  1810. C                                                                      C
  1811. C======================================================================C
  1812.       include 'UIFlow.com'
  1813.       include 'UIFlow.indx'
  1814. C
  1815.       do 11 j   = 2, jmax1
  1816.           ioff  = (j-1) * imaxl + ibeg
  1817.        do 10 i  = 2, imax1
  1818.          ij     = i + ioff
  1819.         amu(ij) = airv * ( t(ij) / airt ) ** 1.5 * ( airsum ) /
  1820.      1                   ( t(ij) + airs )
  1821. 10     continue
  1822. 11    continue
  1823. C
  1824.       return
  1825.       end
  1826. C======================================================================C
  1827.       subroutine lbndry
  1828. C                                                                      C
  1829. C     Purpose:  Prescribe boundary conditions at the zi-minus (left)   C
  1830. C               boundary of the calculation domain.                    C
  1831. C                                                                      C
  1832. C======================================================================C
  1833.       include 'UIFlow.com'
  1834.       igrl = ngrid
  1835.       include 'UIFlow.indx'
  1836. C
  1837.       i = 1
  1838.       do 11 n = 1, nsxm
  1839. C
  1840.        wmxm(n) = wmair
  1841.        if (model .eq. 0) rhxm(n) = rhgs
  1842.        if (model .eq. 1) rhxm(n) = pref * wmair / ( gascon * txm(n) )
  1843.        if (model .eq. 2) then
  1844.         yn2xm   = 1.0 - fuxm(n)  - o2xm(n) - h2oxm(n) - co2xm(n)
  1845.         rwmxm   = fuxm(n)  / wmfu  + o2xm(n)  / wmox + yn2xm / wmn2
  1846.      1          + co2xm(n) / wmco2 + h2oxm(n) / wmh2o
  1847.         wmxm(n) = 1.0 / rwmxm
  1848.         rhxm(n) = pref * wmxm(n) / (gascon * txm(n) + 1.0e-30)
  1849.        endif
  1850. C
  1851.           jfl  = jfxm(n,ngrid)
  1852.       jll  = jlxm(n,ngrid)
  1853.        do 10 j = jfl,jll
  1854.          ijf   = i + (j-1) * imaxl + ibeg
  1855.         rav        = 0.5 * ( r(ijf) + r(ijf-imaxl) )
  1856.         u   (ijf)  = ubxm(n)
  1857.         v   (ijf)  = vbxm(n)
  1858.         w   (ijf)  = wbxm(n)
  1859.         p   (ijf)  = 0.0
  1860.         t   (ijf)  = txm (n)
  1861.         tke (ijf)  = tkxm(n)
  1862.         tde (ijf)  = tdxm(n)
  1863.         f   (ijf)  = fxm(n)
  1864.         g   (ijf)  = gxm(n)
  1865.         yfu (ijf)  = fuxm(n)
  1866.         yo2 (ijf)  = o2xm(n)
  1867.         yh2o(ijf)  = h2oxm(n)
  1868.         yco2(ijf)  = co2xm(n)
  1869.         yn2 (ijf)  = yn2xm
  1870.         gam (ijf)  = vscty
  1871.         amu (ijf)  = vscty
  1872.         amut(ijf)  = cd * rhxm(n) * tkxm(n) * tkxm(n)
  1873.      1               / ( tdxm(n) + 1.0e-30 )
  1874.         wmol(ijf)  = wmxm(n)
  1875.         rho (ijf)  = rhxm(n)
  1876.         cu  (ijf)  = a11(ijf)*u(ijf) + a12(ijf)*v(ijf)
  1877.         cx  (ijf)  = cu (ijf)*rho(ijf)
  1878. 10     continue
  1879. 11    continue
  1880. C
  1881.       return
  1882.       end
  1883. C======================================================================C
  1884.       subroutine masfrc (igrl)
  1885. C                                                                      C
  1886. C     Purpose:  Calculate the mass fractions of species based on law   C
  1887. C               of conservation of atoms.                              C
  1888. C                                                                      C
  1889. C======================================================================C
  1890.       include 'UIFlow.com'
  1891.       include 'UIFlow.indx'
  1892. C
  1893.       rwmfu  = 1.0 / wmfu
  1894.       rwmox  = 1.0 / wmox
  1895.       rwmco2 = 1.0 / wmco2
  1896.       rwmh2o = 1.0 / wmh2o
  1897.       rwmn2  = 1.0 / wmn2
  1898. C
  1899.       do 11 j    = 2, jmax1
  1900.           ioff   = (j-1) * imaxl + ibeg
  1901.        do 10 i   = 2, imax1
  1902.            ij    = i + ioff
  1903.         fmyfu    = f(ij) - yfu(ij)
  1904.         yo2(ij)  = rox * ( 1.0 - f(ij) ) - rat(3) * fmyfu
  1905.         yco2(ij) = rat(1) * fmyfu
  1906.         yh2o(ij) = rat(2) * fmyfu
  1907.         yo2(ij)  = amax1 ( yo2(ij) , 1.0e-20 )
  1908.         yco2(ij) = amax1 ( yco2(ij), 0.0 )
  1909.         yh2o(ij) = amax1 ( yh2o(ij), 0.0 )
  1910.         yn2(ij)  = 1.0 - yfu(ij) - yo2(ij) - yco2(ij) - yh2o(ij)
  1911.         rwmol    = yfu(ij)  * rwmfu   +  yo2(ij)  * rwmox
  1912.      1           + yco2(ij) * rwmco2  +  yh2o(ij) * rwmh2o
  1913.      2           + yn2(ij)  * rwmn2
  1914.         wmol(ij) = 1.0 / rwmol
  1915. 10     continue
  1916. 11    continue
  1917. C
  1918.       return
  1919.       end
  1920. C======================================================================C
  1921.       subroutine moment (igrf)
  1922. C                                                                      C
  1923. C     Purpose:  Perform multi-grid V cycle on continuity and           C
  1924. C               momentum equations.                                    C
  1925. C                                                                                                                   C
  1926. C======================================================================C
  1927.       include 'UIFlow.com'
  1928.       igrl = ngrid
  1929.       include 'UIFlow.indx'
  1930. C
  1931.       flxin = 0.0
  1932. C
  1933.       do 100 n = 1, nsxm
  1934.        if ( kbxm(n) .eq. 2 ) then
  1935.         i = 1
  1936.         do 10 j = jfxm(n,ngrid), jlxm(n,ngrid)
  1937.             ij  = i + ibeg + (j-1) * imaxl
  1938.          flxin  = flxin + cx(ij)
  1939. 10      continue
  1940.        endif
  1941. 100   continue
  1942. C
  1943.       do 200 n = 1, nsxp
  1944.        if ( kbxp(n) .eq. 2 ) then
  1945.         i = imax1
  1946.         do 20 j = jfxp(n,ngrid), jlxp(n,ngrid)
  1947.             ij  = i + ibeg + (j-1) * imaxl
  1948.          flxin  = flxin + cx(ij)
  1949. 20      continue
  1950.        endif
  1951. 200   continue
  1952. C
  1953.       do 300 n = 1, nsym
  1954.        if ( kbym(n) .eq. 2 ) then
  1955.         j = 1
  1956.         do 30 i = ifym(n,ngrid), ilym(n,ngrid)
  1957.             ij  = i + ibeg + (j-1) * imaxl
  1958.          flxin  = flxin + cy(ij)
  1959. 30      continue
  1960.        endif
  1961. 300   continue
  1962. C
  1963.       do 400 n = 1, nsyp
  1964.        if ( kbyp(n) .eq. 2 ) then
  1965.         j = jmax1
  1966.         do 40 i = ifyp(n,ngrid), ilyp(n,ngrid)
  1967.             ij  = i + ibeg + (j-1) * imaxl
  1968.          flxin  = flxin + cy(ij)
  1969. 40      continue
  1970.        endif
  1971. 400   continue
  1972. C
  1973.       return
  1974.       end
  1975. C======================================================================C
  1976.       subroutine mssout(igrl)
  1977. C                                                                      C
  1978. C     Purpose:  Ensure that integral mass flow out of the domain is    C
  1979. C               equal to the influx of mass.                           C
  1980. C                                                                      C
  1981. C======================================================================C
  1982.       include 'UIFlow.com'
  1983.       include 'UIFlow.indx'
  1984. C
  1985.       floact = 1.0e-30
  1986.       flxrto = 0.0
  1987. C
  1988.       i = imax1
  1989.       do 20 j = 2, jmax1
  1990.            ij = i + ibeg + (j-1) * imaxl
  1991.        floact = floact + abs( cx(ij) )
  1992. 20    continue
  1993. C
  1994.       flxrto  = flxin / floact
  1995. C
  1996.       do 21 j = 2, jmax1
  1997.            ij = i + ibeg + (j-1) * imaxl
  1998.        cx(ij) = cx(ij) * flxrto
  1999. 21    continue
  2000. C
  2001.       return
  2002.       end
  2003. C======================================================================C
  2004.       subroutine next (igrf, igrl)
  2005. C                                                                      C
  2006. C     Purpose:  Determine the next step to be performed                C
  2007. C               in the V-cycle.                                        C
  2008. C                                                                      C
  2009. C======================================================================C
  2010.       include 'UIFlow.com'
  2011. C
  2012.       prolng  =  .false.
  2013.       restrc  =  .false.
  2014.       relax   =  .false.
  2015.       niterv  =  nitr(igrl)
  2016. C
  2017.       if (niter .lt. niterv) then
  2018.        relax = .true.
  2019.        return
  2020.       else
  2021.        if(igrf .eq. 1) return
  2022.        if(igrl .eq. 1) upleg = .true.
  2023.        if (upleg) prolng = .true.
  2024.        if (.not. upleg) restrc = .true.
  2025.       endif
  2026. C
  2027.       return
  2028.       end
  2029. C======================================================================C
  2030.       subroutine onestp(igrl)
  2031. C                                                                      C
  2032. C     Purpose:  Calculate all quantities pertinent to the simulation   C
  2033. C               of a turbulent premixed flame.                         C
  2034. C                                                                      C
  2035. C======================================================================C
  2036.       include 'UIFlow.com'
  2037.       include 'UIFlow.indx'
  2038. C
  2039.       call scalar (igrl, 8, f  )
  2040.       call scalar (igrl, 9, yfu)
  2041.       call masfrc (igrl)
  2042.       call trap   (igrl)
  2043.       call scalar (igrl, 5, h  )
  2044.       call props  (igrl)
  2045. C
  2046.       return
  2047.       end
  2048. C======================================================================C
  2049.       subroutine outp
  2050. C                                                                      C
  2051. C     Purpose:  Write out input information to serve as an echo check. C
  2052. C                                                                      C
  2053. C======================================================================C
  2054.       include 'UIFlow.com'
  2055. C
  2056. C              Problem description
  2057. C
  2058.       write (8,1002) klam, kcomp, knorth, kplax, kswrl, kvisc, model, 
  2059.      1               ksc,  kadj
  2060.  1002 format(5x, 'Laminar Flow index (1=Laminar, 0=Turbulent) = ', i5/
  2061.      1       5x, 'Compressible Flow Index (1=Compressible )   = ', i5/
  2062.      2       5x, 'Non-Orthogonal Coordinate index (0 = Orth)  = ', i5/
  2063.      3       5x, 'Planar/Axisymm index (1 = planar)           = ', i5/
  2064.      4       5x, 'Swirl flow index (0 = no swirl)             = ', i5/
  2065.      5       5x, 'Viscosity model(0=fixed laminar viscosity)  = ', i5/
  2066.      6       5x, 'Density model                               = ', i5/
  2067.      7       5x, 'Solve for scalar variable (0 = do not solve)= ', i5/
  2068.      5       5x, 'Block-correction Index(0=no corrections)    = ', i5/)
  2069. C
  2070. C              Grid data
  2071. C
  2072.       write (8,1003) ngrid, ncelx, ncely
  2073.  1003 format(5x, 'Number of Levels of grids                   = ', i5/
  2074.      1       5x, 'Number of Cells (finest grid) in x-dirn     = ', i5/
  2075.      2       5x, 'Number of Cells (finest grid) in y-dirn     = ', i5/)
  2076. C
  2077. C            Segment data
  2078. C
  2079.       write (8,1004) nsxm, nsxp, nsym, nsyp
  2080.  1004 format(5x, 'Number of Segments on x-minus Boundary      = ', i5/
  2081.      1       5x, 'Number of Segments on x-plus  Boundary      = ', i5/
  2082.      2       5x, 'Number of Segments on y-minus Boundary      = ', i5/
  2083.      3       5x, 'Number of Segments on y-plus  Boundary      = ', i5/)
  2084. C
  2085. C            Flow variables Data
  2086. C
  2087.       do 10 n = 1, nsxm
  2088.       write (8,1005) n, jfxm(n,ngrid), jlxm(n,ngrid), kbxm(n)
  2089.       write (8,1006) ubxm(n), vbxm(n), wbxm(n), tkxm(n), tdxm(n), 
  2090.      1               rhxm(n), fxm (n), gxm (n), fuxm(n), txm (n), 
  2091.      2               sxm (n)
  2092.   10  continue
  2093.  1005 format (5x,'Segment Number = ', i5/
  2094.      1        5x,'  jfxm = ',i3,'  jlxm = ',i3/
  2095.      3        5x,'Boundary Index       = ', i5/)
  2096.  1006 format (5x,'u velocity           = ', e12.3/
  2097.      1        5x,'v velocity           = ', e12.3/
  2098.      2        5x,'w velocity           = ', e12.3/
  2099.      3        5x,'turb. kinetic energy = ', e12.3/
  2100.      4        5x,'turb dissipation rate= ', e12.3/
  2101.      5        5x,'density              = ', e12.3/
  2102.      6        5x,'mixture fraction     = ', e12.3/
  2103.      7        5x,'conc. fluctuation    = ', e12.3/
  2104.      8        5x,'fuel fraction        = ', e12.3/
  2105.      9        5x,'temperature          = ', e12.3/
  2106.      9        5x,'general scalar       = ', e12.3/)
  2107. C
  2108.       do 20 n = 1, nsxp
  2109.       write (8,1007) n, jfxp(n,ngrid), jlxp(n,ngrid), kbxp(n)
  2110.       write (8,1006) ubxp(n), vbxp(n), wbxp(n), tkxp(n), tdxp(n), 
  2111.      1               rhxp(n), fxp (n), gxp (n), fuxp(n), txp (n), 
  2112.      2               sxp (n)
  2113.   20  continue
  2114.  1007 format (5x,'Segment Number = ', i5/
  2115.      1        5x,'  jfxp = ',i3,'  jlxp = ',i3/
  2116.      3        5x,'Boundary Index       = ', i5/)
  2117. C
  2118.       do 30 n = 1, nsym
  2119.       write (8,1008) n, ifym(n,ngrid), ilym(n,ngrid), kbym(n)
  2120.       write (8,1006) ubym(n), vbym(n), wbym(n), tkym(n), tdym(n), 
  2121.      1               rhym(n), fym (n), gym (n), fuym(n), tym (n), 
  2122.      2               sym (n)
  2123.   30  continue
  2124.  1008 format (5x,'Segment Number = ', i5/
  2125.      1        5x,'  ifym = ',i3,'  ilym = ',i3/
  2126.      3        5x,'Boundary Index       = ', i5/)
  2127. C
  2128.       do 40 n = 1, nsyp
  2129.       write (8,1009) n, ifyp(n,ngrid), ilyp(n,ngrid), kbyp(n)
  2130.       write (8,1006) ubyp(n), vbyp(n), wbyp(n), tkyp(n), tdyp(n), 
  2131.      1               rhyp(n), fyp (n), gyp (n), fuyp(n), typ (n), 
  2132.      2               syp (n)
  2133.   40  continue
  2134.  1009 format (5x,'Segment Number = ', i5/
  2135.      1        5x,'  ifyp = ',i3,'  ilyp = ',i3/
  2136.      3        5x,'Boundary Index       = ', i5/)
  2137. C
  2138. C                    Guessed Values
  2139. C
  2140.       write (8,1012)
  2141.  1012 format(5x, 'Guessed Values'/5x,15(1h-)/)
  2142.       write (8,1006) ugs, vgs, wgs, tkgs, tdgs, rhgs, fgs, ggs, 
  2143.      1              fugs, tgs,sgs
  2144. C
  2145. C                 Print and Solution Parameters
  2146. C
  2147.       write (8,1013) (iprint(i), i=1,18)
  2148.  1013 format (5x, 'iprint array =', 18i5/)
  2149.       write (8,1014) (relx(n), n=1,11)
  2150.  1014 format (5x, 'relax factor array =', 1p11e8.1/)
  2151.       write (8,1015) (nswp(n), n=1,13)
  2152.  1015 format (5x, 'number of sweeps   =', 13i5/)
  2153.       write (8,1016) (prl(n), n=1,13)
  2154.       write (8,1017) (prt(n), n=1,13)
  2155.  1016 format (5x, 'Laminar Prandtl numbers =', 13f6.1)
  2156.  1017 format (5x, 'Turb.Prandtl numbers    =', 13f6.1)
  2157.       write (8,1018) (tolr(n), n=1,ngrid)
  2158.  1018 format (5x, 'error tolerances', 5e8.3)
  2159. C
  2160. C                    Reference values
  2161. C
  2162.       write(8,1019) refu, refv, refw, refc
  2163.  1019 format (5x, 'Reference values '/
  2164.      1        5x, 'u-velocity             =', e12.3/
  2165.      2        5x, 'v-velocity             =', e12.3/
  2166.      3        5x, 'w-velocity             =', e12.3/     
  2167.      4        5x, 'pressure               =', e12.3/)
  2168.        return
  2169.        end
  2170. C======================================================================C
  2171.       subroutine output(igrl)
  2172. C                                                                      C
  2173. C     Purpose:  Print out flow variables after the calculation has     C
  2174. C               terminated.                                            C
  2175. C                                                                      C
  2176. C======================================================================C
  2177.       include 'UIFlow.com'
  2178. C
  2179.       imaxl  =  imax (igrl)
  2180.       jmaxl  =  jmax (igrl)
  2181.       ibeg   =  nbeg (igrl)
  2182. C
  2183.       if (iprint(1).eq.1) write (8, 1001)
  2184.       if (iprint(1).eq.1) call  plane (imaxl, jmaxl,ibeg, u)
  2185. C
  2186.       if (iprint(2).eq.1) write (8, 1002)
  2187.       if (iprint(2).eq.1) call  plane (imaxl, jmaxl,ibeg, v)
  2188. C
  2189.       if (iprint(3).eq.1) write (8, 1003)
  2190.       if (iprint(3).eq.1) call  plane (imaxl, jmaxl,ibeg, p)
  2191. C
  2192.       if (iprint(4).eq.1) write (8, 1004)
  2193.       if (iprint(4).eq.1) call  plane (imaxl, jmaxl,ibeg, w)
  2194. C
  2195.       if (iprint(5).eq.1) write (8, 1005)
  2196.       if (iprint(5).eq.1) call  plane (imaxl, jmaxl,ibeg, h)
  2197. C
  2198.       if (iprint(5).eq.1) write (8, 1012)
  2199.       if (iprint(5).eq.1) call  plane (imaxl, jmaxl,ibeg, t)
  2200. C
  2201.       if (iprint(6).eq.1) write (8, 1006)
  2202.       if (iprint(6).eq.1) call  plane (imaxl, jmaxl,ibeg, tke)
  2203. C
  2204.       if (iprint(6).eq.1) write (8, 1013)
  2205.       if (iprint(6).eq.1) call  plane (imaxl, jmaxl,ibeg, amut)
  2206. C
  2207.       if (iprint(7).eq.1) write (8, 1007)
  2208.       if (iprint(7).eq.1) call  plane (imaxl, jmaxl,ibeg, tde)
  2209. C
  2210.       if (iprint(8).eq.1) write (8, 1008)
  2211.       if (iprint(8).eq.1) call  plane (imaxl, jmaxl,ibeg, f)
  2212. C
  2213.       if (iprint(9).eq.1) write (8, 1009)
  2214.       if (iprint(9).eq.1) call  plane (imaxl, jmaxl,ibeg, yfu)
  2215. C
  2216.       if (iprint(10).eq.1) write (8, 1010)
  2217.       if (iprint(10).eq.1) call  plane (imaxl, jmaxl,ibeg, g)
  2218. C
  2219.       if (iprint(11).eq.1) write (8, 1011)
  2220.       if (iprint(11).eq.1) call  plane (imaxl, jmaxl,ibeg, rho)
  2221. C
  2222.       if (iprint(12).eq.1) write (8, 1014)
  2223.       if (iprint(12).eq.1) call  plane (imaxl, jmaxl,ibeg, yco2)
  2224. C
  2225.       if (iprint(12).eq.1) write (8, 1015)
  2226.       if (iprint(12).eq.1) call  plane (imaxl, jmaxl,ibeg, yh2o)
  2227. C
  2228.       if (iprint(12).eq.1) write (8, 1016)
  2229.       if (iprint(12).eq.1) call  plane (imaxl, jmaxl,ibeg, yo2)
  2230. C
  2231.       if (iprint(12).eq.1) write (8, 1017)
  2232.       if (iprint(12).eq.1) call  plane (imaxl, jmaxl,ibeg, yn2)
  2233. C
  2234. 1001  format (// 5x, 20(1H-), '    u - Velocity   ', 20(1H-)/)
  2235. 1002  format (// 5x, 20(1H-), '    v - Velocity   ', 20(1H-)/)
  2236. 1003  format (// 5x, 20(1H-), '      Pressure     ', 20(1H-)/)
  2237. 1004  format (// 5x, 20(1H-), '    w - Velocity   ', 20(1H-)/)
  2238. 1005  format (// 5x, 20(1H-), '      Enthalpy     ', 20(1H-)/)
  2239. 1006  format (// 5x, 20(1H-), '    Kinetic Energy ', 20(1H-)/)
  2240. 1007  format (// 5x, 20(1H-), '     Dissipation   ', 20(1H-)/)
  2241. 1008  format (// 5x, 20(1H-), '    Mix.Fraction   ', 20(1H-)/)
  2242. 1009  format (// 5x, 20(1H-), '     Fuel Frac.    ', 20(1H-)/)
  2243. 1010  format (// 5x, 20(1H-), ' Conc. Fluctuation ', 20(1H-)/)
  2244. 1011  format (// 5x, 20(1H-), '      Density      ', 20(1H-)/)
  2245. 1012  format (// 5x, 20(1H-), '     Temperature   ', 20(1H-)/)
  2246. 1013  format (// 5x, 20(1H-), '     Turb. Visc.   ', 20(1H-)/)
  2247. 1014  format (// 5x, 20(1H-), '        CO2        ', 20(1H-)/)
  2248. 1015  format (// 5x, 20(1H-), '        H2O        ', 20(1H-)/)
  2249. 1016  format (// 5x, 20(1H-), '      Oxygen       ', 20(1H-)/)
  2250. 1017  format (// 5x, 20(1H-), '     Nitrogen      ', 20(1H-)/)
  2251. C
  2252.       return
  2253.       end
  2254. C======================================================================C
  2255.       subroutine pcoefc (igrl)
  2256. C                                                                      C
  2257. C     Purpose:  Add terms to the pressure correction equation which    C
  2258. C               account for density fluctuations.  This is only        C
  2259. C               called if the flow is compressible.                    C
  2260. C                                                                      C
  2261. C======================================================================C
  2262.       include 'UIFlow.com'
  2263.       include 'UIFlow.indx'
  2264. C
  2265.       do 11 j  = 2, jmax1
  2266.         ioff   = (j-1) * imaxl + ibeg
  2267.         ioffm  = ioff - imaxl
  2268.         ioffp  = ioff + imaxl
  2269.        do 10 i = 2, imax1
  2270.           ij   = i + ioff
  2271.           ijm  = i + ioffm
  2272.           ijp  = i + ioffp
  2273.           ij1  = ij - 1
  2274.           ijp1 = ij + 1
  2275.         ae(ij) = ae(ij) + amax1( -cu(ij), 0. ) * wmol(ijp1)
  2276.      1                  / ( gamma * gascon  * t(ijp1) )
  2277.         aw(ij) = aw(ij) + amax1( cu(ij1), 0. )* wmol(ij1)
  2278.      1                  / ( gamma * gascon * t(ij-1) )
  2279.         an(ij) = an(ij) + amax1( -cv(ij), 0. )* wmol(ijp)
  2280.      1                  / ( gamma * gascon * t(ijp) )
  2281.         as(ij) = as(ij) + amax1( cv(ijm), 0. )* wmol(ijm)
  2282.      1                  / ( gamma * gascon * t(ijm) )
  2283.         ap(ij) = ap(ij) + ( amax1( cu(ij), 0.) + amax1( -cu(ij1), 0.)
  2284.      1                  +   amax1( cv(ij), 0.) + amax1( -cv(ijm), 0.) )
  2285.      2                  *   wmol(ij)/( gamma * gascon * t(ij) )
  2286. 10     continue
  2287. 11    continue
  2288. C
  2289.       return
  2290.       end
  2291. C======================================================================C
  2292.       subroutine pcoefo (igrl)
  2293. C                                                                      C
  2294. C     Purpose:  Add terms to the pressure correction source term       C
  2295. C               which account for use of a highly NON-ORTHOGONAL grid. C
  2296. C                                                                      C
  2297. C======================================================================C
  2298.       include 'UIFlow.com'
  2299.       include 'UIFlow.indx'
  2300. C
  2301.       call bcor   (igrl, pp)
  2302.       call grad   (igrl, pp)
  2303.       call cpgrad (igrl)
  2304. C
  2305.       do 10 j  = 2, jmax1
  2306.          ioff  = (j-1) * imaxl + ibeg 
  2307.        do 11 i = 2, imax1
  2308.           ij   = i + ioff
  2309.         sp(ij) = x1(ij-1) - x1(ij) + x2(ij-imaxl) - x2(ij)
  2310. 11     continue
  2311. 10    continue
  2312. C
  2313.       return
  2314.       end
  2315. C======================================================================C
  2316.       subroutine plane (idim, jdim, nbg, varl)
  2317. C                                                                      C
  2318. C     Purpose:  Print out a given field variable.                      C
  2319. C                                                                      C
  2320. C======================================================================C
  2321.       dimension varl (*)
  2322.       ip = 1
  2323.       jp = 1
  2324.       istp = ip * 8
  2325. C
  2326.       do 10 ilo = 1,idim,istp
  2327.             ihi = min0 (ilo+istp-1, idim)
  2328.        write(8,1001)
  2329.        do 20 jj = 1,jdim,jp
  2330.               j = jdim - jj + 1
  2331.             ijl = ilo + (j-1) * idim + nbg  
  2332.             ijh = ihi + (j-1) * idim + nbg
  2333.         write(8, 1002) j, (varl(ij), ij = ijl,ijh,ip)
  2334. 20     continue
  2335. 10    continue
  2336. C
  2337. 1001  format(4x, 'j')
  2338. 1002  format (2x, i3, 3x, 1p8e10.3)
  2339.       return
  2340.       end
  2341. C======================================================================C
  2342.       subroutine prlong (igrf)
  2343. C                                                                      C
  2344. C     Purpose:  Call subroutines to prolongate variables from one      C
  2345. C               fine grid to the next.  Used in the nested iteration   C
  2346. C               phase of the FMG cycle.                                C
  2347. C                                                                      C
  2348. C======================================================================C
  2349.       include 'UIFlow.com'
  2350. C
  2351.       call  prolcx (igrf)
  2352.       call  prolcy (igrf)
  2353.       call  prolcu (igrf)
  2354.       call  prolcv (igrf)
  2355.       call  prols  (igrf, u, 0)
  2356.       call  prols  (igrf, v, 0)
  2357.       call  prols  (igrf, p, 1)
  2358. C
  2359.       if (klam .eq. 0) then
  2360.        call  prols (igrf, gam,  0)
  2361.        call  prols (igrf, tke,  0)
  2362.        call  prols (igrf, tde,  0)
  2363.        call  prols (igrf, amut, 0)
  2364.        if (model .eq. 3) then
  2365.         call prols (igrf, g, 0)
  2366.        endif
  2367.       endif
  2368. C
  2369.       if (kswrl.eq. 1) call  prols (igrf, w, 0)
  2370. C
  2371.       if (model .gt. 0) call  prols (igrf, h,    0)
  2372.       if (model .gt. 0) call  prols (igrf, t,    0)
  2373.       if (model .gt. 0) call  prols (igrf, rho,  0)
  2374.       if (model .gt. 1) call  prols (igrf, f,    0)
  2375.       if (model .gt. 1) call  prols (igrf, yfu,  0)
  2376.       if (model .gt. 1) call  prols (igrf, yo2,  0)
  2377.       if (model .gt. 1) call  prols (igrf, yco2, 0)
  2378.       if (model .gt. 1) call  prols (igrf, wmol, 0)
  2379.       if (model .eq. 2) call  prols (igrf, yn2,  0)
  2380.       if (model .eq. 2) call  prols (igrf, yh2o, 0)
  2381. C
  2382.       return
  2383.       end
  2384. C======================================================================C
  2385.       subroutine prodn (igrl)
  2386. C                                                                      C
  2387. C     Purpose:  Calculate the production of turbulence kinetic energy. C
  2388. C                                                                      C
  2389. C======================================================================C
  2390.       include 'UIFlow.com'
  2391.       include 'UIFlow.indx'
  2392. C
  2393.       do 11 j  = 2, jmax1
  2394.          ioff  = (j-1) * imaxl + ibeg
  2395.         ioffm  = ioff - imaxl
  2396.         ioffp  = ioff + imaxl
  2397.        do 10 i = 2, imax1
  2398.             ij = i + ioff
  2399.            imj = ij - 1
  2400.            ipj = ij + 1
  2401.            ijm = i + ioffm
  2402.            ijp = i + ioffp
  2403. C
  2404. C.... Compute derivatives in the zi and eta system.
  2405. C
  2406.         dudzi  = fx(ij) * u(ipj) + fx1(ij) * u(ij) - fx(imj) * u(ij)
  2407.      1         - fx1(imj) * u(imj)
  2408.         dvdzi  = fx(ij) * v(ipj) + fx1(ij) * v(ij) - fx(imj) * v(ij)
  2409.      1         - fx1(imj) * v(imj)
  2410.         dudet  = fy(ij) * u(ijp) + fy1(ij) * u(ij) - fy(ijm)
  2411.      1         * u(ij) - fy1(ijm) * u(ijm)
  2412.         dvdet  = fy(ij) * v(ijp) + fy1(ij) * v(ij) - fy(ijm)
  2413.      1         * v(ij) - fy1(ijm) * v(ijm)
  2414. C
  2415. C.... Compute derivatives in the x and y system.
  2416. C
  2417.         rajb   = 1.0 / ajb(ij)
  2418. C
  2419.         dudx   = ( dudzi * dux(ij) + dudet * duy(ij) ) * rajb
  2420.         dvdx   = ( dvdzi * dux(ij) + dvdet * duy(ij) ) * rajb
  2421.         dudy   = ( dudzi * dvx(ij) + dudet * dvy(ij) ) * rajb
  2422.         dvdy   = ( dvdzi * dvx(ij) + dvdet * dvy(ij) ) * rajb
  2423. C              
  2424.         prod(ij) = 2.0 * ( dudx * dudx  +  dvdy * dvdy )
  2425.      1           + (dudy + dvdx)**2
  2426. 10     continue
  2427. 11    continue
  2428. C
  2429.       if ( kplax .eq. 0 ) then
  2430.        do 21 j    = 2, jmax1
  2431.           ioff    = (j-1) * imaxl + ibeg
  2432.          ioffm    = ioff - imaxl
  2433.         do 20 i   = 2, imax1
  2434.              ij   = i + ioff
  2435.             ij1   = i + ioffm
  2436.          rav      = 0.25 * ( r(ij) + r(ij-1) + r(ij1) + r(ij1-1) )
  2437.          prod(ij) = prod(ij) + 2.0 * ( v(ij) / rav )**2
  2438. 20      continue
  2439. 21     continue
  2440.       endif
  2441. C
  2442.       do 31 j    = 2, jmax1
  2443.          ioff    = (j-1) * imaxl + ibeg
  2444.        do 30 i   = 2, imax1
  2445.             ij   = i + ioff
  2446.         prod(ij) = prod(ij) * amut(ij)
  2447. 30     continue
  2448. 31    continue
  2449. C
  2450. C.... Modify production term if calculating a swirl flow.
  2451. C
  2452.       if ( kswrl .eq. 1 ) then
  2453.        do 41 j  = 2, jmax1
  2454.           ioff  = (j-1) * imaxl + ibeg
  2455.          ioffm  = ioff - imaxl
  2456.         do 40 i = 2, imax1
  2457.              ij = i + ioff
  2458.             imj = ij - 1
  2459.             ij1 = i + ioffm
  2460.          rav    = 0.25 * ( r(ij) + r(imj) + r(ij1) + r(ij1-1) )
  2461.          rajb   = 1.0 / ajb(ij)
  2462.          dwdzi  = fx(ij) * w(ij+1) + fx1(ij) * w(ij) - fx(imj) * w(ij)
  2463.      1          - fx1(imj) * w(imj)
  2464.          dwdet  = fy(ij) * w(ij+imaxl) + fy1(ij) * w(ij) - fy(ij1)
  2465.      1            * w(ij) - fy1(ij1) * w(ij1)
  2466.          dwdx   = ( dwdzi * dux(ij) + dwdet * duy(ij) ) * rajb
  2467.          dwdy   = ( dwdzi * dvx(ij) + dwdet * dvy(ij) ) * rajb
  2468.          prod(ij) = prod(ij) + amut(ij) * ( dwdx * dwdx
  2469.      1            + ( dwdy - w(ij)/rav )**2 )
  2470. 40      continue
  2471. 41     continue
  2472.       endif
  2473. C
  2474.       return
  2475.       end
  2476. C======================================================================C
  2477.       subroutine prolcu (igrl)
  2478. C                                                                      C
  2479. C     Purpose:  Prolongate corrections to contravariant velocities     C
  2480. C               perpendicular to lines of constant zi.                 C
  2481. C                                                                      C
  2482. C======================================================================C
  2483.       include 'UIFlow.com'
  2484. C
  2485.       igrlp  =  igrl + 1
  2486.       imaxc  =  imax (igrl)
  2487.       jmaxc  =  jmax (igrl)
  2488.       imaxf  =  imax (igrlp)
  2489.       jmaxf  =  jmax (igrlp)
  2490.       imaxc1 =  imaxc - 1
  2491.       jmaxc1 =  jmaxc - 1
  2492.       ibegc  =  nbeg (igrl)
  2493.       ibegf  =  nbeg (igrlp)
  2494. C
  2495. C.... Initialise corrections.
  2496. C
  2497.       do 11 jc  = 1, jmaxc 
  2498.        do 10 ic = 1, imaxc 
  2499.            ijc  = ic + (jc - 1) * imaxc + ibegc
  2500.         c(ijc)  = 0.0
  2501. 10     continue
  2502. 11    continue
  2503. C
  2504. C.... Evaluate corrections on the coarse grid.
  2505. C
  2506.       do 21 jc  = 2, jmaxc1
  2507.        do 20 ic = 1, imaxc1
  2508.            ijc  = ic + (jc - 1) * imaxc + ibegc
  2509.        ijf  = iru (ijc)
  2510.         c(ijc)  = cu(ijc) - ( cu(ijf) + cu(ijf-imaxf) )
  2511.         c(ijc)  = 0.5 * c(ijc)
  2512. 20     continue
  2513. 21    continue
  2514. C
  2515. C.... Prolongate corrections to the finer grid.
  2516. C
  2517.       a1 = 0.75
  2518.       a2 = 0.25
  2519. C
  2520.       do 31 ic  = 2, imaxc1
  2521.        do 30 jc = 2, jmaxc1
  2522.             ijc = ic + (jc - 1) * imaxc + ibegc
  2523.         ijf = iru (ijc)
  2524.            ijfm = ijf - imaxf
  2525.        ijc1 = ijc - 1
  2526.            ijcp = ijc + imaxc
  2527.            ijcm = ijc - imaxc
  2528.         cu(ijf)    = cu(ijf)  + a1 * c(ijc) + a2 * c(ijcp)
  2529.         cu(ijfm)   = cu(ijfm) + a1 * c(ijc) + a2 * c(ijcm)
  2530.         cu(ijf-1)  = cu(ijf-1)  + 0.5 * ( a1*( c(ijc)+c(ijc1) )
  2531.      1             + a2 * ( c(ijcp)+c(ijc1+imaxc) ) )
  2532.         cu(ijfm-1) = cu(ijfm-1) + 0.5 * ( a1*( c(ijc)+c(ijc1) )
  2533.      1             + a2 * ( c(ijcm)+c(ijc1-imaxc) ) )
  2534. 30     continue
  2535. 31    continue
  2536. C
  2537.       return
  2538.       end
  2539. C======================================================================C
  2540.       subroutine prolcv (igrl)
  2541. C                                                                      C
  2542. C     Purpose:  Prolongate corrections to contravariant velocities     C
  2543. C               perpendicular to lines of constant eta.                C
  2544. C                                                                      C
  2545. C======================================================================C
  2546.       include 'UIFlow.com'
  2547. C
  2548.       igrlp  =  igrl + 1
  2549.       imaxc  =  imax (igrl)
  2550.       jmaxc  =  jmax (igrl)
  2551.       imaxf  =  imax (igrlp)
  2552.       jmaxf  =  jmax (igrlp)
  2553.       imaxc1 =  imaxc - 1
  2554.       jmaxc1 =  jmaxc - 1
  2555.       ibegc  =  nbeg (igrl)
  2556.       ibegf  =  nbeg (igrlp)
  2557. C
  2558. C.... Initialise corrections.
  2559. C
  2560.       do 10 jc  = 1, jmaxc 
  2561.        do 11 ic = 1, imaxc 
  2562.            ijc  = ic + (jc - 1) * imaxc + ibegc
  2563.         c(ijc)  = 0.0
  2564. 11     continue
  2565. 10    continue
  2566. C
  2567. C.... Evaluate corrections on the coarse grid.
  2568. C
  2569.       do 20 jc  = 1, jmaxc1
  2570.        do 21 ic = 2, imaxc1
  2571.            ijc  = ic + (jc - 1) * imaxc + ibegc
  2572.        ijf  = iru (ijc)
  2573.         c(ijc)  = cv(ijc) - ( cv(ijf) + cv(ijf-1) )
  2574.         c(ijc)  = 0.5 * c(ijc)
  2575. 21     continue
  2576. 20    continue
  2577. C
  2578. C.... Prolongate corrections to the fine grid.
  2579. C
  2580.       a1 = 0.75
  2581.       a2 = 0.25
  2582. C
  2583.       do 30 ic  = 2, imaxc1
  2584.        do 31 jc = 2, jmaxc1
  2585.             ijc = ic + (jc - 1) * imaxc + ibegc
  2586.        ijc1 = ijc - imaxc
  2587.            ijcm = ijc - 1
  2588.            ijcp = ijc + 1
  2589.         ijf = iru (ijc)
  2590.            ijf1 = ijf - 1
  2591.            ijfm = ijf - imaxf
  2592.         cv(ijf)  = cv(ijf)  + a1 * c(ijc) + a2 * c(ijcp)
  2593.         cv(ijf1) = cv(ijf1) + a1 * c(ijc) + a2 * c(ijcm)
  2594.         cv(ijfm) = cv(ijfm) + 0.5 * ( a1*( c(ijc)+c(ijc1) )
  2595.      1           + a2*( c(ijcp)+c(ijc1+1) ) )
  2596.         cv(ijfm-1) = cv(ijfm-1) + 0.5 * ( a1*( c(ijc)+c(ijc1) )
  2597.      1             + a2*( c(ijcm)+c(ijc1-1) ) )
  2598. 31     continue
  2599. 30    continue
  2600. C
  2601.       return
  2602.       end
  2603. C======================================================================C
  2604.       subroutine prolcx (igrl)
  2605. C                                                                      C
  2606. C     Purpose:  Prolongate corrections to mass fluxes perpendicular    C
  2607. C               to lines of constant zi.                               C
  2608. C                                                                      C
  2609. C======================================================================C
  2610.       include 'UIFlow.com'
  2611. C
  2612.       igrlp  =  igrl + 1
  2613.       imaxc  =  imax (igrl)
  2614.       jmaxc  =  jmax (igrl)
  2615.       imaxf  =  imax (igrlp)
  2616.       jmaxf  =  jmax (igrlp)
  2617.       imaxc1 =  imaxc - 1
  2618.       jmaxc1 =  jmaxc - 1
  2619.       ibegc  =  nbeg (igrl)
  2620.       ibegf  =  nbeg (igrlp)
  2621. C
  2622. C.... Initialise corrections.
  2623. C
  2624.       do 11 jc  = 1, jmaxc 
  2625.        do 10 ic = 1, imaxc 
  2626.            ijc  = ic + (jc - 1) * imaxc + ibegc
  2627.         c(ijc)  = 0.0
  2628. 10     continue
  2629. 11    continue
  2630. C
  2631. C.... Evaluate corrections on the coarse grid.
  2632. C
  2633.       do 21 jc  = 2, jmaxc1
  2634.        do 20 ic = 1, imaxc1
  2635.            ijc  = ic + (jc - 1) * imaxc + ibegc
  2636.        ijf  = iru (ijc)
  2637.         c(ijc)  = cx(ijc) - ( cx(ijf) + cx(ijf-imaxf) )
  2638.         c(ijc)  = 0.5 * c(ijc)
  2639. 20     continue
  2640. 21    continue
  2641. C
  2642. C.... Prolongate corrections to the finer grid.
  2643. C
  2644.       a1 = 0.75
  2645.       a2 = 0.25
  2646. C
  2647.       do 31 ic  = 2, imaxc1
  2648.        do 30 jc = 2, jmaxc1
  2649.             ijc = ic + (jc - 1) * imaxc + ibegc
  2650.         ijf = iru (ijc)
  2651.            ijfm = ijf - imaxf
  2652.        ijc1 = ijc - 1
  2653.            ijcp = ijc + imaxc
  2654.            ijcm = ijc - imaxc
  2655.         cx(ijf)    = cx(ijf)  + a1 * c(ijc) + a2 * c(ijcp)
  2656.         cx(ijfm)   = cx(ijfm) + a1 * c(ijc) + a2 * c(ijcm)
  2657.         cx(ijf-1)  = cx(ijf-1)  + 0.5 * ( a1*( c(ijc)+c(ijc1) )
  2658.      1             + a2 * ( c(ijcp)+c(ijc1+imaxc) ) )
  2659.         cx(ijfm-1) = cx(ijfm-1) + 0.5 * ( a1*( c(ijc)+c(ijc1) )
  2660.      1             + a2 * ( c(ijcm)+c(ijc1-imaxc) ) )
  2661. 30     continue
  2662. 31    continue
  2663. C
  2664.       return
  2665.       end
  2666. C======================================================================C
  2667.       subroutine prolcy (igrl)
  2668. C                                                                      C
  2669. C     Purpose:  Prolongate corrections to mass fluxes perpendicular    C
  2670. C               to lines of constant eta.                              C
  2671. C                                                                      C
  2672. C======================================================================C
  2673.       include 'UIFlow.com'
  2674. C
  2675.       igrlp  =  igrl + 1
  2676.       imaxc  =  imax (igrl)
  2677.       jmaxc  =  jmax (igrl)
  2678.       imaxf  =  imax (igrlp)
  2679.       jmaxf  =  jmax (igrlp)
  2680.       imaxc1 =  imaxc - 1
  2681.       jmaxc1 =  jmaxc - 1
  2682.       ibegc  =  nbeg (igrl)
  2683.       ibegf  =  nbeg (igrlp)
  2684. C
  2685. C.... Initialise corrections.
  2686. C
  2687.       do 10 jc  = 1, jmaxc 
  2688.        do 11 ic = 1, imaxc 
  2689.            ijc  = ic + (jc - 1) * imaxc + ibegc
  2690.         c(ijc)  = 0.0
  2691. 11     continue
  2692. 10    continue
  2693. C
  2694. C.... Evaluate corrections on the coarse grid.
  2695. C
  2696.       do 20 jc  = 1, jmaxc1
  2697.        do 21 ic = 2, imaxc1
  2698.            ijc  = ic + (jc - 1) * imaxc + ibegc
  2699.        ijf  = iru (ijc)
  2700.         c(ijc)  = cy(ijc) - ( cy(ijf) + cy(ijf-1) )
  2701.         c(ijc)  = 0.5 * c(ijc)
  2702. 21     continue
  2703. 20    continue
  2704. C
  2705. C.... Prolongate corrections to the fine grid.
  2706. C
  2707.       a1 = 0.75
  2708.       a2 = 0.25
  2709. C
  2710.       do 30 ic  = 2, imaxc1
  2711.        do 31 jc = 2, jmaxc1
  2712.             ijc = ic + (jc - 1) * imaxc + ibegc
  2713.        ijc1 = ijc - imaxc
  2714.            ijcm = ijc - 1
  2715.            ijcp = ijc + 1
  2716.         ijf = iru (ijc)
  2717.            ijf1 = ijf - 1
  2718.            ijfm = ijf - imaxf
  2719.         cy(ijf)  = cy(ijf)  + a1 * c(ijc) + a2 * c(ijcp)
  2720.         cy(ijf1) = cy(ijf1) + a1 * c(ijc) + a2 * c(ijcm)
  2721.         cy(ijfm) = cy(ijfm) + 0.5 * ( a1*( c(ijc)+c(ijc1) )
  2722.      1           + a2*( c(ijcp)+c(ijc1+1) ) )
  2723.         cy(ijfm-1) = cy(ijfm-1) + 0.5 * ( a1*( c(ijc)+c(ijc1) )
  2724.      1             + a2*( c(ijcm)+c(ijc1-1) ) )
  2725. 31     continue
  2726. 30    continue
  2727. C
  2728.       return
  2729.       end
  2730. C======================================================================C
  2731.       subroutine prols (igrl,q1,indx)
  2732. C                                                                      C
  2733. C     Purpose:  Prolongate corrections to variable (q1) located        C
  2734. C               at cell center.                                        C
  2735. C                                                                      C
  2736. C======================================================================C
  2737.       include 'UIFlow.com'
  2738.       dimension q1(*)
  2739. C
  2740.       igrlp  =  igrl + 1
  2741.       imaxc  =  imax (igrl)
  2742.       jmaxc  =  jmax (igrl)
  2743.       imaxf  =  imax (igrlp)
  2744.       jmaxf  =  jmax (igrlp)
  2745.       imaxc1 =  imaxc - 1
  2746.       jmaxc1 =  jmaxc - 1
  2747.       ibegc  =  nbeg (igrl)
  2748.       ibegf  =  nbeg (igrlp)
  2749. C
  2750. C.... Initialise corrections.
  2751. C
  2752.       do 10 jc  = 1, jmaxc
  2753.        do 11 ic = 1, imaxc
  2754.           ijc   = ic + (jc-1) *imaxc + ibegc
  2755.         c(ijc)  = 0.0
  2756. 11     continue
  2757. 10    continue
  2758. C
  2759. C.... Evaluate corrections on the coarse grid.
  2760. C
  2761.       do 20 jc = 2, jmaxc1
  2762.        do 21 ic = 2, imaxc1
  2763.            ijc = ic + (jc - 1) * imaxc + ibegc
  2764.        ijf = iru (ijc)
  2765.         c(ijc) = q1(ijc) - 0.25 * ( q1(ijf)       + q1(ijf-1)
  2766.      1                            + q1(ijf-imaxf) + q1(ijf-imaxf-1) )
  2767. 21     continue
  2768. 20    continue
  2769. C
  2770. C.... Extrapolate corrections at boundaries.
  2771. C
  2772.       if(indx .eq. 1) call bcor (igrl,c)
  2773. C
  2774. C.... Prolongate corrections to the fine grid.
  2775. C
  2776.       a1 = 9./16.
  2777.       a2 = 3./16.
  2778.       a3 = 3./16.
  2779.       a4 = 1./16.
  2780. C
  2781.       do 30 ic  = 2, imaxc1
  2782.        do 31 jc = 2, jmaxc1
  2783.            ijc  = ic + (jc - 1) * imaxc + ibegc
  2784.           ijcm  = ijc - imaxc
  2785.           ijcp  = ijc + imaxc
  2786.        ijf  = iru (ijc)
  2787.           ijf1  = ijf - 1
  2788.           ijfm  = ijf - imaxf
  2789.         q1(ijf)    = q1(ijf) + a1 * c(ijc)  + a2 * c(ijc+1)
  2790.      1                       + a3 * c(ijcp) + a4 * c(ijcp+1)
  2791.         q1(ijf1)   = q1(ijf1) + a1 * c(ijc)  + a2 * c(ijc-1)
  2792.      1                        + a3 * c(ijcp) + a4 * c(ijcp-1)
  2793.         q1(ijfm)   = q1(ijfm) + a1 * c(ijc)  + a2 * c(ijc+1)
  2794.      1                        + a3 * c(ijcm) + a4 * c(ijcm+1)
  2795.         q1(ijfm-1) = q1(ijfm-1) + a1 * c(ijc)  + a2 * c(ijc-1)
  2796.      1                          + a3 * c(ijcm) + a4 * c(ijcm-1)
  2797. 31     continue
  2798. 30    continue
  2799. C
  2800.       return
  2801.       end
  2802. C======================================================================C
  2803.       subroutine props (igrl)
  2804. C
  2805. C     Purpose:  Call subroutines which will calculate the temperature  C
  2806. C               and density for the compressible flow and premixed     C
  2807. C               flame models.                                          C
  2808. C                                                                      C
  2809. C======================================================================C
  2810.       include 'UIFlow.com'
  2811.       include 'UIFlow.indx'
  2812. C
  2813.       call cpenth (igrl)
  2814.       call temp   (igrl)
  2815.       call dens   (igrl)
  2816. C
  2817.       return
  2818.       end
  2819. C======================================================================C
  2820.       subroutine rbndry
  2821. C                                                                      C
  2822. C     Purpose:  Prescribe boundary conditions for the zi plus (right)  C
  2823. C               boundary of the calculation domain.                    C
  2824. C                                                                      C
  2825. C======================================================================C
  2826.       include 'UIFlow.com'
  2827.       igrl = ngrid
  2828.       include 'UIFlow.indx'
  2829. C
  2830.       i = imaxl
  2831. C
  2832.       do 11 n = 1, nsxp
  2833. C
  2834.        wmxp(n) = wmair
  2835.        if (model .eq. 0) rhxp(n) = rhgs
  2836.        if (model .eq. 1) rhxp(n) = pref * wmair / (gascon * txp(n))
  2837.        if (model .eq. 2) then
  2838.         yn2xp   = 1.0 - fuxp(n)  - o2xp(n) - h2oxp(n) - co2xp(n)
  2839.         rwmxp   = fuxp(n)  / wmfu  + o2xp(n)  / wmox + yn2xp / wmn2
  2840.      1          + co2xp(n) / wmco2 + h2oxp(n) / wmh2o
  2841.         wmxp(n) = 1.0 / rwmxp
  2842.         rhxp(n) = pref * wmxp(n) / (gascon * txp(n) + 1.0e-30)
  2843.        endif
  2844. C
  2845.            jfl = jfxp(n,ngrid)
  2846.        jll = jlxp(n,ngrid)
  2847.        do 10 j = jfl,jll
  2848.          ijf   = i + (j-1) * imaxl + ibeg
  2849.          ijf1  = ijf - 1
  2850.         u   (ijf)  = ubxp(n)
  2851.         v   (ijf)  = vbxp(n)
  2852.         w   (ijf)  = wbxp(n)
  2853.         p   (ijf)  = 0.0
  2854.         t   (ijf)  = txp (n)
  2855.         tke (ijf)  = tkxp(n)
  2856.         tde (ijf)  = tdxp(n)
  2857.         f   (ijf)  = fxp(n)
  2858.         g   (ijf)  = gxp(n)
  2859.         yfu (ijf)  = fuxp(n)
  2860.         yo2 (ijf)  = o2xp(n)
  2861.         yh2o(ijf)  = h2oxp(n)
  2862.         yco2(ijf)  = co2xp(n)
  2863.         yn2 (ijf)  = yn2xp
  2864.         gam (ijf)  = vscty
  2865.         amu (ijf)  = vscty
  2866.         amut(ijf)  = cd * rhxp(n) * tkxp(n) * tkxp(n)
  2867.      1               / ( tdxp(n) + 1.0e-30 )
  2868.         wmol(ijf)  = wmxp(n)
  2869.         rho (ijf)  = rhxp(n)
  2870.         cu  (ijf1) = a11(ijf1) * ubxp(n) + a12(ijf1) * vbxp(n)
  2871.         cx  (ijf1) = cu (ijf1) * rho(ijf)
  2872. 10     continue
  2873. 11    continue
  2874. C
  2875.       return
  2876.       end
  2877. C======================================================================C
  2878.       subroutine resid(igrl, q, sum)
  2879. C                                                                      C
  2880. C     Purpose:  Compute the residuals for any given variable.          C
  2881. C                                                                      C
  2882. C======================================================================C
  2883.       include 'UIFlow.com'
  2884.       dimension q(*)
  2885.       include 'UIFlow.indx'
  2886. C
  2887.       sum = 0.0
  2888. C
  2889.       do 11 j   = 2, jmax1
  2890.        do 10 i  = 2, imax1
  2891.            ijf  = i + (j-1) * imaxl + ibeg
  2892.         rs(ijf) = su(ijf) + aw(ijf) * q(ijf-1) + ae(ijf) * q(ijf+1) +
  2893.      1            as(ijf) * q(ijf-imaxl) + an(ijf) * q(ijf+imaxl) -
  2894.      2            ap(ijf) * q(ijf)
  2895.         sum     = sum + abs (rs(ijf))
  2896. 10     continue
  2897. 11    continue
  2898. C
  2899.       return
  2900.       end
  2901. C======================================================================C
  2902.       subroutine restbs (igrl, q)
  2903. C                                                                      C
  2904. C     Purpose:  Restrict boundary values of the given variable.        C
  2905. C                                                                      C
  2906. C======================================================================C
  2907.       include 'UIFlow.com'
  2908.       dimension q(*)
  2909. C
  2910.       igrm   = igrl - 1
  2911.       imaxf  = imax (igrl)
  2912.       jmaxf  = jmax (igrl)
  2913.       ibegf  = nbeg (igrl)
  2914.       imaxc  = imax (igrm)
  2915.       jmaxc  = jmax (igrm)
  2916.       imaxc1 = imaxc - 1
  2917.       jmaxc1 = jmaxc - 1
  2918.       ibegc  = nbeg (igrm)
  2919. C
  2920. C.... Restrict x - minus boundary.
  2921. C
  2922.       i = 1
  2923.       do 10 j = 2, jmaxc1
  2924.           ijc = i + (j-1) * imaxc + ibegc
  2925.           ijf = iru(ijc)
  2926.        q(ijc) = 0.5 * ( q(ijf) + q(ijf-imaxf) )
  2927. 10    continue
  2928. C
  2929. C.... Restrict x - plus boundary.
  2930. C
  2931.       i = imaxc
  2932.       do 11 j = 2, jmaxc1
  2933.           ijc = i + (j-1) * imaxc + ibegc
  2934.           ijf = iru(ijc)
  2935.        q(ijc) = 0.5 * ( q(ijf) + q(ijf-imaxf) )
  2936. 11    continue
  2937. C
  2938. C.... Restrict y - minus boundary.
  2939. C
  2940.       j = 1
  2941.       do 20 i = 2, imaxc1
  2942.           ijc = i + (j-1) * imaxc + ibegc
  2943.           ijf = iru(ijc)
  2944.        q(ijc) = 0.5 * ( q(ijf) + q(ijf-1) )
  2945. 20    continue
  2946. C
  2947. C.... Restrict y - plus boundary.
  2948. C
  2949.       j  =  jmaxc
  2950.       do 21 i = 2, imaxc1
  2951.           ijc = i + (j-1) * imaxc + ibegc
  2952.           ijf = iru(ijc)
  2953.        q(ijc) = 0.5 * ( q(ijf) + q(ijf-1) )
  2954. 21    continue
  2955. C
  2956.       return
  2957.       end
  2958. C======================================================================C
  2959.       subroutine restfl (igrl)
  2960. C                                                                      C
  2961. C     Purpose:  Restrict the volume and mass fluxes.                   C
  2962. C                                                                      C
  2963. C======================================================================C
  2964.       include 'UIFlow.com'
  2965. C
  2966.       igrlm  = igrl - 1
  2967.       imaxf  = imax (igrl)
  2968.       jmaxf  = jmax (igrl)
  2969.       ibegf  = nbeg (igrl)
  2970.       imaxc  = imax (igrlm)
  2971.       jmaxc  = jmax (igrlm)
  2972.       imaxc1 = imaxc - 1
  2973.       jmaxc1 = jmaxc - 1
  2974.       ibegc  = nbeg (igrlm)
  2975. C
  2976. C.... Restrict fluxes perpendicular to lines of constant zi.
  2977. C
  2978.       do 11 j   = 2, jmaxc1
  2979.        do 10 i  = 1, imaxc1
  2980.            ijc  = i + (j-1) * imaxc + ibegc
  2981.            ijf  = iru(ijc)
  2982.         cx(ijc) = cx(ijf) + cx(ijf-imaxf)
  2983.         cu(ijc) = cx(ijc) / (amax1 (sign (rho(ijc)  ,  cx(ijc)), 0.0)
  2984.      1                    +  amax1 (sign (rho(ijc+1), -cx(ijc)), 0.0))
  2985. 10     continue
  2986. 11    continue
  2987. C
  2988. C.... Restrict fluxes perpendicular to lines of constant eta.
  2989. C
  2990.       do 21 j   = 1, jmaxc1
  2991.        do 20 i  = 2, imaxc1
  2992.            ijc  = i + (j-1) * imaxc + ibegc
  2993.            ijf  = iru(ijc)
  2994.         cy(ijc) = cy(ijf) + cy(ijf-1)
  2995.         cv(ijc) = cy(ijc) / (amax1 (sign (rho(ijc),  cy(ijc)), 0.0)
  2996.      1          +  amax1 (sign (rho(ijc+imaxc), -cy(ijc)), 0.0))
  2997. 20     continue
  2998. 21    continue
  2999. C
  3000.       return
  3001.       end
  3002. C======================================================================C
  3003.       subroutine restr (igrl)
  3004. C                                                                      C
  3005. C     Purpose:  Restrict the residuals to a coarser grid.              C
  3006. C                                                                      C
  3007. C======================================================================C
  3008.       include 'UIFlow.com'
  3009. C
  3010.       igrm   =  igrl - 1
  3011.       imaxc  =  imax (igrm)
  3012.       jmaxc  =  jmax (igrm)
  3013.       imaxf  =  imax (igrl)
  3014.       jmaxf  =  jmax (igrl)
  3015.       imaxc1 =  imaxc - 1
  3016.       jmaxc1 =  jmaxc - 1
  3017.       ibegc  =  nbeg (igrm)
  3018.       ibegf  =  nbeg (igrl)
  3019. C
  3020. C.... Calculate diffusion terms and coefficients.
  3021. C
  3022.       call  dflux (igrl)
  3023.       call  dflux (igrm)
  3024.       call  coeff (igrl)
  3025.       call  coeff (igrm)
  3026. C
  3027. C.... Calculate source terms for u equation.
  3028. C               
  3029.       call  trsrc  (igrl, u)
  3030.       call  trsrc  (igrm, u)
  3031.       call  srcu   (igrl)
  3032.       call  srcu   (igrm)
  3033.       call  urelax (igrl, 1,u)
  3034.       call  urelax (igrm, 1,u)
  3035. C
  3036. C.... Calculate residuals and restrict.
  3037. C
  3038.       call  resid (igrl,u,sum)
  3039.       call  resid (igrm,u,sum)
  3040. C
  3041.       do 11 j  = 2, jmaxc1
  3042.        do 10 i = 2, imaxc1
  3043.           ijc  = i + (j-1)*imaxc + ibegc
  3044.       ijf  = iru(ijc)
  3045.          ijfm  = ijf - imaxf
  3046.         resux(ijc) = rs(ijf) + rs(ijf-1) + rs(ijfm) + rs(ijfm-1)
  3047.         resu (ijc) = resux(ijc) - ( rs(ijc) - resu(ijc) )
  3048. 10     continue
  3049. 11    continue
  3050. C
  3051. C.... Calculate source terms for v equation.
  3052. C
  3053.       call  trsrc  (igrl, v)
  3054.       call  trsrc  (igrm, v)
  3055.       call  srcv   (igrl)
  3056.       call  srcv   (igrm)
  3057.       call  urelax (igrl,2,v)
  3058.       call  urelax (igrm,2,v)
  3059. C
  3060. C.... Calculate residuals and restrict.
  3061. C
  3062.       call  resid (igrl,v,sum)
  3063.       call  resid (igrm,v,sum)
  3064. C
  3065.       do 21 j  = 2, jmaxc1
  3066.        do 20 i = 2, imaxc1
  3067.           ijc  = i + (j-1)*imaxc + ibegc
  3068.       ijf  = iru(ijc)
  3069.          ijfm  = ijf - imaxf
  3070.         resvx(ijc) = rs(ijf) + rs(ijf-1) + rs(ijfm) + rs(ijfm-1)
  3071.         resv (ijc) = resvx(ijc) - ( rs(ijc) - resv(ijc) )
  3072. 20     continue
  3073. 21    continue
  3074. C
  3075.       return
  3076.       end
  3077. C======================================================================C
  3078.       subroutine rests (igrl,q)
  3079. C                                                                      C
  3080. C     Purpose:  Restrict cell centered quantity.                       C
  3081. C                                                                      C
  3082. C======================================================================C
  3083.       include 'UIFlow.com'
  3084.       dimension q(*)
  3085. C
  3086.       igrm   = igrl - 1
  3087.       imaxf  = imax (igrl)
  3088.       jmaxf  = jmax (igrl)
  3089.       ibegf  = nbeg (igrl)
  3090.       imaxc  = imax (igrm)
  3091.       jmaxc  = jmax (igrm)
  3092.       imaxc1 = imaxc - 1
  3093.       jmaxc1 = jmaxc - 1
  3094.       ibegc  = nbeg (igrm)
  3095. C
  3096. C.... Restrict internal variables.
  3097. C
  3098.       do 11 j  = 2, jmaxc1 
  3099.        do 10 i = 2, imaxc1 
  3100.           ijc  = i + (j-1) * imaxc + ibegc
  3101.           ijf1 = iru(ijc)
  3102.       ijf2 = ijf1 - 1
  3103.       ijf3 = ijf1 - imaxf
  3104.       ijf4 = ijf1 - imaxf  - 1
  3105.         q(ijc) = 0.25 * ( q(ijf1) + q(ijf2) + q(ijf3) + q(ijf4) )
  3106. 10     continue
  3107. 11    continue
  3108. C
  3109. C.... Restrict boundary values.
  3110. C
  3111.       call restbs (igrl, q)
  3112. C
  3113.       return
  3114.       end
  3115. C======================================================================C
  3116.       subroutine restv (igrl)
  3117. C                                                                      C
  3118. C     Purpose:  Restrict values of cell centered quantities as used    C
  3119. C               in the V cycle.                                        C
  3120. C                                                                      C
  3121. C======================================================================C
  3122.       include 'UIFlow.com'
  3123. C
  3124.       call rests (igrl, u)
  3125.       call rests (igrl, v)
  3126.       call rests (igrl, p)
  3127.       call rests (igrl, amu)
  3128.       call rests (igrl, gam)
  3129.       call rests (igrl, rho)
  3130. C
  3131.       if (kswrl .eq. 1) call rests(igrl, w)
  3132.       if (klam  .eq. 0) call rests(igrl, amut)
  3133. C
  3134.       return
  3135.       end
  3136. C======================================================================C
  3137.       subroutine scalar (igrf, nv, q)
  3138. C                                                                      C
  3139. C     Purpose:  Calculate coefficients and obtain solution of the      C
  3140. C               given variable.                                        C
  3141. C                                                                      C
  3142. C======================================================================C
  3143.       include 'UIFlow.com'
  3144.       dimension q(*)
  3145.       igrl  = igrf
  3146.       include 'UIFlow.indx'
  3147. C
  3148.       call visc  (igrl, nv)
  3149.       call dflux (igrl)
  3150.       call coeff (igrl)
  3151.       call trsrc (igrl, q)
  3152.       call sorce (igrl, nv)
  3153.       call urelax(igrl, nv, q)
  3154.       call solve (igrl, nv, q)
  3155.       call resid (igrl, q, sum)
  3156.       error (igrl, nv) = sum
  3157. C
  3158.       return
  3159.       end
  3160. C======================================================================C
  3161.       subroutine scalrs (igrf)
  3162. C                                                                      C
  3163. C     Purpose:  Solve appropriate scalars depending on the problem     C
  3164. C               description.                                           C
  3165. C                                                                      C
  3166. C======================================================================C
  3167.       include 'UIFlow.com'
  3168.       igrl = igrf
  3169.       include 'UIFlow.indx'
  3170. C
  3171.       do 10 itnke = 1, nitke
  3172. C
  3173.        if (kswrl .eq. 1) call scalar (igrf, 4, w  )
  3174. C
  3175.        if (klam .eq. 0) then
  3176.         call prodn  (igrf)
  3177.         call scalar (igrf, 6, tke)
  3178.         call scalar (igrf, 7, tde)
  3179.        endif
  3180. C
  3181.        if (model .eq. 1) call scalar (igrf, 5, h)
  3182.        if (model .eq. 2) call onestp (igrf)
  3183.        if (model .eq. 3) call fstchm (igrf)
  3184. 10    continue
  3185. C
  3186.       return
  3187.       end
  3188. C======================================================================C
  3189.       subroutine search
  3190. C                                                                      C
  3191. C     Purpose:  Calculate oxygen mass fractions at all inlets based on C
  3192. C               a prescribed fuel mass fraction.  This precludes the   C
  3193. C               specification of an already burned fuel rich mixture   C
  3194. C               at the inlet.                                          C
  3195. C                                                                      C
  3196. C======================================================================C
  3197.       include 'UIFlow.com'
  3198. C
  3199.       do 10 n = 1, nsxm
  3200.        if ( kbxm(n) .eq. 2 ) then
  3201.         yair     = 1.0 - fuxm(n)
  3202.         o2xm (n) = rox * yair
  3203.         h2oxm(n) = 0.0
  3204.         co2xm(n) = 0.0
  3205.        endif
  3206. 10    continue
  3207. C
  3208.       do 20 n = 1, nsxp
  3209.        if ( kbxp(n) .eq. 2 ) then
  3210.         yair     = 1.0 - fuxp(n)
  3211.         o2xp (n) = rox * yair
  3212.         h2oxp(n) = 0.0
  3213.         co2xp(n) = 0.0
  3214.        endif
  3215. 20    continue
  3216. C
  3217.       do 30 n = 1, nsym
  3218.        if ( kbym(n) .eq. 2 ) then
  3219.         yair     = 1.0 - fuym(n)
  3220.         o2ym (n) = rox * yair
  3221.         h2oym(n) = 0.0
  3222.         co2ym(n) = 0.0
  3223.        endif
  3224. 30    continue
  3225. C
  3226.       do 40 n = 1, nsyp
  3227.        if ( kbyp(n) .eq. 2 ) then
  3228.         yair     = 1.0 - fuyp(n)
  3229.         o2yp (n) = rox * yair
  3230.         h2oyp(n) = 0.0
  3231.         co2yp(n) = 0.0
  3232.        endif
  3233. 40    continue
  3234. C
  3235.       return
  3236.       end
  3237. C======================================================================C
  3238.       subroutine solve (igrl, nv, phi)
  3239. C                                                                      C
  3240. C     Purpose:  Obtain solution to the given flow variable by use of   C
  3241. C               an ADI technique.                                      C
  3242. C                                                                      C
  3243. C======================================================================C
  3244.       include 'UIFlow.com'
  3245.       dimension phi(*), ad(200), bd(200), cc(200), dd(200)
  3246.       include 'UIFlow.indx'
  3247. C
  3248.       ifrst = 2
  3249.       ilst  = imax1
  3250.       jfrst = 2
  3251.       jlst  = jmax1
  3252.       iswps = 0
  3253. C
  3254. 100   continue
  3255. C
  3256. C.... Set up tri-diagonal matrix for an EAST-WEST inversion.
  3257. C
  3258.       do 21 j  = jfrst, jlst
  3259.        do 20 i = ifrst, ilst
  3260.             ij = ibeg + i + (j-1) * imaxl
  3261.            ijp = ij + imaxl
  3262.            ijm = ij - imaxl
  3263.         bd(i) = - aw(ij)
  3264.         dd(i) =   ap(ij) - beta * an(ij)
  3265.         ad(i) = - ae(ij)
  3266.         cc(i) =   an(ij) * ( phi(ijp) - beta * phi(ij) ) +
  3267.      1            as(ij) *   phi(ijm) + su(ij) + sp(ij)
  3268. 20     continue
  3269.         ijf = ifrst + ibeg + (j-1) * imaxl
  3270.         ijl = ilst  + ibeg + (j-1) * imaxl
  3271.         cc(ifrst) = cc(ifrst) + aw(ijf) * phi(ijf-1)
  3272.         cc(ilst)  = cc(ilst)  + ae(ijl) * phi(ijl+1)
  3273. C
  3274. C.... Perform forward elimination.
  3275. C
  3276.       do 22 k = ifrst + 1, ilst
  3277.            k1 = k - 1
  3278.        dd(k)  = dd(k) - ( bd(k) * ad(k1) ) / dd(k1)
  3279.        cc(k)  = cc(k) - ( bd(k) * cc(k1) ) / dd(k1)
  3280. 22    continue
  3281. C
  3282. C.... Back substitute.
  3283. C
  3284.       ij = ilst  + ibeg + (j-1) * imaxl
  3285.       phi(ij) = cc(ilst) / dd(ilst)
  3286. C
  3287.       do 23 k  = 1, ilst - 2
  3288.             kk = ilst - k
  3289.            ij  = kk + ibeg + (j-1) * imaxl
  3290.        phi(ij) = ( cc(kk) - ad(kk) * phi(ij+1) ) / dd(kk)
  3291. 23    continue
  3292. C
  3293. 21    continue
  3294. C
  3295. C.... Call pcoefo if the grid structure is highly non-orthogonal.
  3296. C
  3297.       if ( (nv .eq. 3) .and. (knorth .eq. 1) ) call pcoefo (igrl)
  3298. C
  3299. C.... Set up tri-diagonal matrix for a NORTH-SOUTH inversion.
  3300. C
  3301.       do 31 i  = ifrst, ilst
  3302.        do 30 j = jfrst, jlst
  3303.             ij = ibeg + i + (j-1) * imaxl
  3304.         bd(j)  = - as(ij)
  3305.         dd(j)  =   ap(ij) - beta * ae(ij)
  3306.         ad(j)  = - an(ij)
  3307.         cc(j)  =   ae(ij) * ( phi(ij+1) - beta * phi(ij) ) +
  3308.      1             aw(ij) *   phi(ij-1) + su(ij) + sp(ij)
  3309. 30     continue
  3310.         ijf = i + ibeg + (jfrst-1) * imaxl
  3311.         ij1 = ijf - imaxl
  3312.         ijl = i + ibeg + (jlst-1) * imaxl
  3313.         ij2 = ijl + imaxl
  3314.         cc(jfrst) = cc(jfrst) + as(ijf) * phi(ij1)
  3315.         cc(jlst)  = cc(jlst)  + an(ijl) * phi(ij2)
  3316. C
  3317. C.... Perform forward elimination.
  3318. C
  3319.       do 32 k = jfrst + 1, jlst
  3320.            k1 = k - 1
  3321.        dd(k)  = dd(k) - ( bd(k) * ad(k1) ) / dd(k1)
  3322.        cc(k)  = cc(k) - ( bd(k) * cc(k1) ) / dd(k1)
  3323. 32    continue
  3324. C
  3325. C.... Back substitute.
  3326. C
  3327.       ij = i + ibeg + (jlst-1) * imaxl
  3328.       phi(ij) = cc(jlst) / dd(jlst)
  3329. C
  3330.       do 33 k  = 1, jlst - 2
  3331.            kk  = jlst - k
  3332.           ij   = i + ibeg + (kk-1) * imaxl
  3333.           ijp  = ij + imaxl
  3334.        phi(ij) = ( cc(kk) - ad(kk)*phi(ijp) ) / dd(kk)
  3335. 33    continue
  3336. 31    continue
  3337. C
  3338. C.... Call pcoefo if the grid structure is highly non-orthogonal.
  3339. C
  3340.       if ( (nv .eq. 3) .and. (knorth .eq. 1) ) call pcoefo (igrl)
  3341. C
  3342.       iswps = iswps + 1
  3343. C
  3344. C.... Check for convergence of algebraic solver.
  3345. C
  3346.       if (iswps .ge. nswp(nv) ) then
  3347.        goto 200
  3348.       else
  3349.        goto 100
  3350.       endif
  3351. C
  3352. 200   continue
  3353.       return
  3354.       end
  3355. C======================================================================C
  3356.       subroutine sorce (igrl, nv)
  3357. C                                                                      C
  3358. C     Purpose:  Call appropriate routine for the calculation of the    C
  3359. C               given variable's source term.                          C
  3360. C                                                                      C
  3361. C======================================================================C
  3362.       include 'UIFlow.com'
  3363.       include 'UIFlow.indx'
  3364. C
  3365.       goto (10, 20, 30, 40, 50, 60, 70, 80, 90, 100), nv
  3366. C
  3367. C.... u - momentum equation, nv = 1
  3368. C
  3369. 10    continue
  3370.       call srcu (igrl)
  3371.       return
  3372. C
  3373. C.... v - momentum equation, nv = 2
  3374. C
  3375. 20    continue
  3376.       call srcv (igrl)
  3377.       return
  3378. C
  3379. C.... p prime equation , nv = 3
  3380. C
  3381. 30    continue
  3382.       return
  3383. C
  3384. C.... w - momentum equation, nv = 4
  3385. C
  3386. 40    continue
  3387.       call srcw (igrl)
  3388.       return
  3389. C
  3390. C.... energy equation, nv = 5
  3391. C
  3392. 50    continue
  3393.       call srch (igrl)
  3394.       return
  3395. C
  3396. C.... turbulence kinetic energy, nv = 6
  3397. C
  3398. 60    continue
  3399.       call  srck (igrl)
  3400.       return
  3401. C
  3402. C.... turbulence dissipation rate, nv = 7
  3403. C
  3404. 70    continue
  3405.       call srcd (igrl)
  3406.       return
  3407. C
  3408. C.... mixture fraction ( f ) , nv = 8
  3409. C
  3410. 80    continue
  3411.       call srcf (igrl)
  3412.       return
  3413. C
  3414. C.... fuel mass fraction ( yfu ), nv = 9
  3415. C
  3416. 90    continue
  3417.       call srcfu (igrl)
  3418.       return
  3419. C
  3420. C.... concentraction fluctuation ( g ), nv = 10
  3421. C
  3422. 100   continue
  3423.       call srcg (igrl)
  3424.       return
  3425.       end
  3426. C======================================================================C
  3427.       subroutine srcd (igrl)
  3428. C                                                                      C
  3429. C     Purpose:  Compute the source term for the dissipation equation.  C
  3430. C               srcd also implements the wall functions by calling     C
  3431. C               walld for each wall segment.                           C
  3432. C                                                                      C
  3433. C======================================================================C
  3434.       include 'UIFlow.com'
  3435.       include 'UIFlow.indx'
  3436. C
  3437. C.... Compute production and dissipation.
  3438. C
  3439.       do 11 j  = 2, jmax1
  3440.          ioff  = (j-1) * imaxl + ibeg
  3441.        do 10 i = 2, imax1
  3442.           ij   = i + ioff
  3443.         sp(ij) = - ce2 * rho (ij) * tde(ij) / tke(ij) * ajb(ij)
  3444.      1         +   amin1( su(ij)/tde(ij), 0.0 )
  3445.         su(ij) =   ce1 * prod(ij) * tde(ij) / tke(ij) * ajb(ij)
  3446.      1         +   amax1( su(ij), 0.0 )
  3447. 10     continue
  3448. 11    continue
  3449. C
  3450. C.... Enforce wall functions on turb. dissipation rate.
  3451. C
  3452.       i = 2
  3453.       do 20 n = 1, nsxm
  3454.        if ( kbxm(n) .ne. 1 ) goto 20
  3455.           jfl  = jfxm(n,igrl)
  3456.       jll  = jlxm(n,igrl)
  3457.        call  walld (igrl,i,i,jfl,jll,xxic)
  3458. 20    continue
  3459. C
  3460.       i = imax1
  3461.       do 30 n = 1, nsxp
  3462.        if ( kbxp(n) .ne. 1 ) goto 30
  3463.            jfl = jfxp(n,igrl)
  3464.        jll = jlxp(n,igrl)
  3465.        call  walld (igrl,i,i,jfl,jll,xxic)
  3466. 30    continue
  3467. C
  3468.       j = 2
  3469.       do 40 n = 1, nsym
  3470.        if ( kbym(n) .ne. 1 ) goto 40
  3471.           ifl  = ifym(n,igrl)
  3472.       ill  = ilym(n,igrl)
  3473.        call  walld (igrl,ifl,ill,j,j,yetac)
  3474. 40    continue
  3475. C
  3476.       j = jmax1
  3477.       do 50 n = 1, nsyp
  3478.        if ( kbyp(n) .ne. 1 ) goto 50
  3479.           ifl  = ifyp(n,igrl)
  3480.       ill  = ilyp(n,igrl)
  3481.        call  walld (igrl,ifl,ill,j,j,yetac)
  3482. 50    continue
  3483. C
  3484.       return
  3485.       end
  3486. C======================================================================C
  3487.       subroutine srcf (igrl)
  3488. C                                                                      C
  3489. C     Purpose:  Compute the source term for the mixture fraction       C
  3490. C               equation.                                              C
  3491. C                                                                      C
  3492. C======================================================================C
  3493.       include 'UIFlow.com'
  3494.       include 'UIFlow.indx'
  3495. C
  3496.       do 11 j  = 2, jmax1
  3497.          ioff  = (j-1) * imaxl + ibeg
  3498.        do 10 i = 2, imax1
  3499.            ij  = i + ioff
  3500.         sp(ij) = amin1( su(ij)/f(ij), 0.0 )
  3501.         su(ij) = amax1( su(ij), 0.0 )
  3502. 10     continue
  3503. 11    continue
  3504. C
  3505.       return
  3506.       end
  3507. C======================================================================C
  3508.       subroutine srcfu (igrl)
  3509. C                                                                      C
  3510. C     Purpose:  Compute the source terms for the fuel fraction         C
  3511. C               equation.                                              C
  3512. C                                                                      C
  3513. C======================================================================C
  3514.       include 'UIFlow.com'
  3515.       include 'UIFlow.indx'
  3516. C
  3517.       if ( klam .eq. 1 ) then
  3518. C
  3519. C.... Source term determination for a LAMINAR flow.
  3520. C
  3521.        do 11 j   = 2, jmax1
  3522.           ioff   = (j-1) * imaxl + ibeg
  3523.         do 10 i  = 2, imax1
  3524.              ij  = i + ioff
  3525.          sp(ij)  = amin1( su(ij)/yfu(ij), 0.0 )
  3526.          su(ij)  = amax1( su(ij), 0.0 )
  3527.          arhn    = - prexp(1) * exp(-acten(1) / t(ij))
  3528.          sor     = arhn * ( yfu(ij)**aa(1) ) * ( yo2(ij)**bb(1) )
  3529.      1           * ( rho(ij)**ab(1) )
  3530.          sp (ij) = sp(ij) + sor * ajb(ij) / yfu(ij)
  3531.          sfu(ij) = sor * ajb(ij)
  3532. 10      continue
  3533. 11     continue
  3534. C
  3535.       elseif ( klam .eq. 0 ) then
  3536. C
  3537. C.... Source term determination for a TURBULENT flow.
  3538. C
  3539.        do 21 j   = 2, jmax1
  3540.           ioff   = (j-1) * imaxl + ibeg
  3541.         do 20 i  = 2, imax1
  3542.              ij  = i + ioff
  3543.          sp(ij)  = amin1( su(ij)/yfu(ij), 0.0 )
  3544.          su(ij)  = amax1( su(ij), 0.0 )
  3545.          arhn    = - prexp(1) * exp( -acten(1) / t(ij) )
  3546.          sor1    = arhn * ( yfu(ij)**aa(1) ) * ( yo2(ij)**bb(1) )
  3547.      1           * ( rho(ij) ** ab(1) )
  3548.          sor2    = - cr * yfu(ij) * rho(ij) * tde(ij) / tke(ij)
  3549.          sor     = amax1 ( sor1, sor2 ) * ajb(ij)
  3550.          sp(ij)  = sp(ij) + sor / yfu(ij)
  3551.          sfu(ij) = sor
  3552. 20      continue
  3553. 21     continue
  3554.       endif
  3555. C
  3556.       return
  3557.       end
  3558. C======================================================================C
  3559.       subroutine srcg (igrl)
  3560. C                                                                      C
  3561. C     Purpose:  Compute the source terms for the concentration         C
  3562. C               fluctuation equation.                                  C
  3563. C                                                                      C
  3564. C======================================================================C
  3565.       include 'UIFlow.com'
  3566.       include 'UIFlow.indx'
  3567. C
  3568.       do 11 j  = 2, jmax1
  3569.          ioff  = (j-1) * imaxl + ibeg
  3570.        do 10 i = 2, imax1
  3571.             ij = i + ioff
  3572.            ij1 = ij - 1
  3573.            ijp = ij + imaxl
  3574.            ijm = ij - imaxl
  3575. C
  3576. C.... Compute derivatives in computational domain.
  3577. C
  3578.         dfdzi  = fx(ij)  * f(ij+1) + fx1(ij)  * f(ij)
  3579.      1         - fx(ij1) * f(ij)   - fx1(ij1) * f(ij1)
  3580.         dfdeta = fy(ij)  * f(ijp)  + fy1(ij)  * f(ij)
  3581.      1         - fy(ijm) * f(ij)   - fy1(ijm) * f(ijm)
  3582. C
  3583. C.... Compute derivatives in the physical space.
  3584. C
  3585.         rajb = 1.0 / ajb(ij)
  3586.         dfdx = ( dfdzi * dux(ij) + dfdeta * duy(ij) ) * rajb
  3587.         dfdy = ( dfdzi * dvx(ij) + dfdeta * dvy(ij) ) * rajb
  3588. C
  3589. C.... Calculate source term for the scalar fluctuation equation.
  3590. C
  3591.         sp(ij) = - cg2 * rho(ij) * tde(ij) * ajb(ij) / tke(ij)
  3592.      1           + amin1( su(ij)/g(ij) , 0.0 )
  3593.         su(ij) = cg2 * amut(ij) * ajb(ij) / prt(10)
  3594.      1         * ( dfdx * dfdx + dfdy * dfdy )
  3595.      2         + amax1( su(ij) , 0.0 )
  3596. 10     continue
  3597. 11    continue
  3598. C
  3599.       return
  3600.       end
  3601. C======================================================================C
  3602.       subroutine srch (igrl)
  3603. C                                                                      C
  3604. C     Purpose:  Calculate source term for the enthalpy equation.       C
  3605. C                                                                      C
  3606. C======================================================================C
  3607.       include 'UIFlow.com'
  3608.       include 'UIFlow.indx'
  3609. C
  3610.       if ( model .ne. 2 ) return
  3611. C
  3612. C.... Source terms which account for specie heat flux.
  3613. C
  3614.       do 11 j  = 2, jmax1
  3615.          ioff  = (j-1) * imaxl + ibeg
  3616.        do 10 i = 2, imax1
  3617.             ij = i + ioff
  3618.         sp(ij) = amin1( su(ij)/h(ij), 0.0 )
  3619.         su(ij) = amax1( su(ij), 0.0 )
  3620.         su(ij) = su(ij) + sfu(ij) * ( rat(4)*hfco2 + rat(5)*hfh2o )
  3621.         sp(ij) = sp(ij) - sfu(ij) * hffu / h(ij)
  3622. 10     continue
  3623. 11    continue
  3624. C
  3625.       return
  3626.       end
  3627. C======================================================================C
  3628.       subroutine srck (igrl)
  3629. C                                                                      C
  3630. C     Purpose:  Compute the source terms in the turbulent kinetic      C
  3631. C               energy equation.                                       C
  3632. C                                                                      C
  3633. C======================================================================C
  3634.       include 'UIFlow.com'
  3635.       include 'UIFlow.indx'
  3636. C
  3637.       do 11 j  = 2, jmax1
  3638.          ioff  = (j-1) * imaxl + ibeg
  3639.        do 10 i = 2, imax1
  3640.           ij   = i + ioff
  3641.         sp(ij) = - rho(ij) * tde(ij) / tke(ij) * ajb(ij)
  3642.      1         +   amin1( su(ij)/tke(ij), 0.0 )
  3643.         su(ij) =   prod(ij) * ajb(ij) + amax1( su(ij), 0.0 )
  3644. 10     continue
  3645. 11    continue
  3646. C
  3647. C.... Modify coefficients and source terms at walls.
  3648. C
  3649.       i = 2
  3650.       do 20 n = 1, nsxm
  3651.        if ( kbxm(n) .eq. 1 ) then
  3652.           jfl = jfxm(n,igrl)
  3653.       jll = jlxm(n,igrl)
  3654.         call  wallk (igrl,i,i,jfl,jll,u,v,w,0.0,0.0,0.0,xxic,aw)
  3655.        endif
  3656. 20    continue
  3657. C
  3658.       i = imax1
  3659.       do 30 n = 1, nsxp
  3660.        if ( kbxp(n) .eq. 1 ) then
  3661.           jfl = jfxp(n,igrl)
  3662.       jll = jlxp(n,igrl)
  3663.         call  wallk (igrl,i,i,jfl,jll,u,v,w,0.0,0.0,0.0,xxic,ae)
  3664.        endif
  3665. 30    continue
  3666. C
  3667.       j = 2
  3668.       do 40 n = 1, nsym
  3669.        if ( kbym(n) .eq. 1 ) then
  3670.           ifl = ifym(n,igrl)
  3671.       ill = ilym(n,igrl)
  3672.         call  wallk (igrl,ifl,ill,j,j,u,v,w,0.0,0.0,0.0,yetac,as)
  3673.        endif
  3674. 40    continue
  3675. C
  3676.       j = jmax1
  3677.       do 50 n = 1, nsyp
  3678.        if ( kbyp(n) .eq. 1 ) then
  3679.           ifl = ifyp(n,igrl)
  3680.       ill = ilyp(n,igrl)
  3681.         call  wallk (igrl,ifl,ill,j,j,u,v,w,0.0,0.0,0.0,yetac,an)
  3682.        endif
  3683. 50    continue
  3684. C
  3685.       return
  3686.       end
  3687. C======================================================================C
  3688.       subroutine srcu (igrl)
  3689. C                                                                      C
  3690. C     Purpose:  Calculate source terms for the u momentum equation.    C
  3691. C                                                                      C
  3692. C======================================================================C
  3693.       include 'UIFlow.com'
  3694.       include 'UIFlow.indx'
  3695. C
  3696. C.... Calculate pressure gradient terms.
  3697. C
  3698.       call grad (igrl, p)
  3699.       call bcor (igrl, dpdx)
  3700.       call bcor (igrl, dpdy)
  3701. C
  3702. C.... Calculate apu and source terms.
  3703. C
  3704.       do 11  j  = 2, jmax1
  3705.           ioff  = (j-1) * imaxl + ibeg
  3706.        do 10 i  = 2, imax1
  3707.             ij  = i + ioff
  3708.         apu(ij) = su(ij)
  3709.         su(ij)  = (dux(ij) * dpdx(ij) + duy(ij) * dpdy(ij)) +
  3710.      1             resu(ij) + su(ij)
  3711.         sp(ij)  = 0.0
  3712. 10     continue
  3713. 11    continue
  3714. C
  3715.       return
  3716.       end
  3717. C======================================================================C
  3718.       subroutine srcv (igrl)
  3719. C                                                                      C
  3720. C     Purpose:  Calculate source terms for the v momentum equation.    C
  3721. C                                                                      C
  3722. C======================================================================C
  3723.       include 'UIFlow.com'
  3724.       include 'UIFlow.indx'
  3725. C
  3726. C.... Calculate source terms axisymmetric flows.
  3727. C
  3728.       if ( kplax .eq. 0 ) then
  3729.        do 11 j  = 2, jmax1
  3730.            ioff = (j-1) * imaxl + ibeg
  3731.         do 10 i = 2, imax1
  3732.             ij  = i + ioff
  3733.            ij1  = ij - imaxl
  3734.          rav    = 0.25 * ( r(ij) + r(ij-1) + r(ij1) + r(ij1-1) )
  3735.          su(ij) = su(ij) + rho(ij) * w(ij) * w(ij) / rav * ajb(ij)
  3736.          sp(ij) = - 2.0 * gam(ij) / (rav * rav) * ajb(ij)
  3737. 10      continue
  3738. 11     continue
  3739.       endif
  3740. C
  3741. C.... Calculate apv and source terms.
  3742. C
  3743.       do 21 j   = 2, jmax1
  3744.           ioff  = (j-1) * imaxl + ibeg
  3745.        do 20 i  = 2, imax1
  3746.             ij  = i + ioff
  3747.         apv(ij) = su(ij)
  3748.         su (ij) = (dvx(ij) * dpdx(ij) + dvy(ij) * dpdy(ij)) +
  3749.      1             resv(ij) + su(ij)
  3750. 20     continue
  3751. 21    continue
  3752. C
  3753.       return
  3754.       end
  3755. C======================================================================C
  3756.       subroutine srcw (igrl)
  3757. C                                                                      C
  3758. C     Purpose:  Calculate source terms for the w momentum equation.    C
  3759. C                                                                      C
  3760. C======================================================================C
  3761.       include 'UIFlow.com'
  3762.       include 'UIFlow.indx'
  3763. C
  3764.       do 11 j  = 2, jmax1
  3765.          ioff  = (j-1) * imaxl + ibeg
  3766.        do 10 i = 2, imax1
  3767.             ij = i  + ioff
  3768.            ij1 = ij - 1
  3769.            ijp = ij + imaxl
  3770.            ijm = ij - imaxl
  3771.         rav    = 0.25 * ( r(ij) + r(ij1) + r(ijm) + r(ijm-1) )
  3772.         dgdzi  = fx(ij)  * gam(ij+1) + fx1(ij)  * gam(ij)
  3773.      1         - fx(ij1) * gam(ij)   - fx1(ij1) * gam(ij1)
  3774.         dgdeta = fy(ij)  * gam(ijp) + fy1(ij)  * gam(ij)
  3775.      1         - fy(ijm) * gam(ij)  - fy1(ijm) * gam(ijm)
  3776.         dgamdr = dgdzi * dvx(ij) + dgdeta * dvy(ij)
  3777.         term1  = - rho(ij) * v(ij) * w(ij) / rav * ajb(ij)
  3778.         term2  = - gam(ij) * w(ij) / ( rav * rav ) * ajb(ij)
  3779.         term3  = - w(ij) * dgamdr  / rav
  3780. C
  3781.         sp(ij) = amin1( su(ij)/w(ij), 0.0 ) + amin1( term1/w(ij), 0.0 )
  3782.      1         + amin1( term2/w(ij) , 0.0 ) + amin1( term3/w(ij), 0.0 )
  3783.         su(ij) = amax1( su(ij), 0.0 ) + amax1( term1, 0.0 )
  3784.      1         + amax1( term2 , 0.0 ) + amax1( term3, 0.0 )
  3785. C
  3786. 10     continue
  3787. 11    continue
  3788.       return
  3789.       end
  3790. C======================================================================C
  3791.       subroutine tbndry
  3792. C                                                                      C
  3793. C     Purpose:  Prescribe boundary conditions on the eta plus (top)    C 
  3794. C               boundary.                                              C
  3795. C                                                                      C
  3796. C======================================================================C
  3797.       include 'UIFlow.com'
  3798.       igrl = ngrid
  3799.       include 'UIFlow.indx'
  3800. C
  3801.       j = jmaxl
  3802. C
  3803.       do 11 n = 1, nsyp
  3804. C
  3805.        wmyp(n) = wmair
  3806.        if (model .eq. 0) rhyp(n) = rhgs
  3807.        if (model .eq. 1) rhyp(n) = pref * wmair / ( gascon * typ(n) )
  3808.        if (model .eq. 2) then
  3809.         yn2yp   = 1.0 - fuyp(n) - o2yp(n) - h2oyp(n) - co2yp(n)
  3810.         rwmyp   = fuyp(n)  / wmfu  + o2yp(n)  / wmox + yn2yp / wmn2
  3811.      1          + co2yp(n) / wmco2 + h2oyp(n) / wmh2o
  3812.         wmyp(n) = 1.0 / rwmyp
  3813.         rhyp(n) = pref * wmyp(n) / (gascon * typ(n) + 1.0e-30)
  3814.        endif
  3815.           ifl = ifyp(n,ngrid)
  3816.       ill = ilyp(n,ngrid)
  3817.        do 10 i = ifl,ill
  3818.          ijf   = i + (j-1) * imaxl + ibeg
  3819.      ijf1  = ijf - imaxl
  3820.         u   (ijf)  = ubyp(n)
  3821.         v   (ijf)  = vbyp(n)
  3822.         w   (ijf)  = wbyp(n)
  3823.         p   (ijf)  = 0.0
  3824.         t   (ijf)  = typ (n)
  3825.         tke (ijf)  = tkyp(n)
  3826.         tde (ijf)  = tdyp(n)
  3827.         f   (ijf)  = fyp(n)
  3828.         g   (ijf)  = gyp(n)
  3829.         yfu (ijf)  = fuyp(n)
  3830.         yo2 (ijf)  = o2yp(n)
  3831.         yh2o(ijf)  = h2oyp(n)
  3832.         yco2(ijf)  = co2yp(n)
  3833.         yn2 (ijf)  = yn2yp
  3834.         gam (ijf)  = vscty
  3835.         amu (ijf)  = vscty
  3836.         amut(ijf)  = cd * rhyp(n) * tkyp(n) * tkyp(n)
  3837.      1               / ( tdyp(n) + 1.0e-30 )
  3838.         wmol(ijf)  = wmyp(n)
  3839.         rho (ijf)  = rhyp(n)
  3840.         cv (ijf1)  = a21(ijf1)*u(ijf) + a22(ijf1)*v(ijf)
  3841.         cy (ijf1)  = cv (ijf1)*rho(ijf)
  3842. 10     continue
  3843. 11    continue
  3844. C
  3845.       return
  3846.       end
  3847. C======================================================================C
  3848.       subroutine temp (igrl)
  3849. C                                                                      C
  3850. C     Purpose:  Calculate fluid temperature for compressible and       C
  3851. C               premixed flame models.                                 C
  3852. C                                                                      C
  3853. C======================================================================C
  3854.       include 'UIFlow.com'
  3855.       include 'UIFlow.indx'
  3856. C
  3857.       relxm = 1.0 - relx(5)
  3858. C
  3859.       if ( model .eq. 1 ) then
  3860. C
  3861. C... Calculate temperature for compressible air flow problem.
  3862. C
  3863.        hmix = rox * hox  +  (1.0 - rox) * hn2
  3864. C
  3865.        do 11 j  = 2, jmax1
  3866.           ioff  = (j-1) * imaxl + ibeg
  3867.         do 10 i = 2, imax1
  3868.              ij = i + ioff
  3869.          t1     = ( h(ij) + hmix ) / cph(ij)
  3870.          t(ij)  = relx(5) * t1 + relxm * t(ij)
  3871. 10      continue
  3872. 11     continue
  3873. C
  3874.       elseif ( model .eq. 2 ) then
  3875. C
  3876.        do 21 j  = 2, jmax1
  3877.           ioff  = (j-1) * imaxl + ibeg
  3878.         do 20 i = 2, imax1
  3879.              ij = i + ioff
  3880.          hmix   =  yfu(ij)  * hfu  +  yo2(ij)  * hox
  3881.      1          +  yn2(ij)  * hn2  +  yco2(ij) * hco2
  3882.      2          +  yh2o(ij) * hh2o
  3883.          t1     = ( h(ij) + hmix ) / cph(ij)
  3884.          t(ij)  = relx(5) * t1  +  relxm * t(ij)
  3885.          t(ij)  = amin1( t(ij), 3500.0 )
  3886. 20      continue
  3887. 21     continue
  3888.       endif 
  3889. C
  3890.       return
  3891.       end
  3892. C======================================================================C
  3893.       subroutine trap (igrl)
  3894. C                                                                      C
  3895. C     Purpose:  Ensure that the sum of all mass fractions is equal     C
  3896. C               to unity.                                              C
  3897. C                                                                      C
  3898. C======================================================================C
  3899.       include 'UIFlow.com'
  3900.       include 'UIFlow.indx'
  3901. C
  3902.       do 11 j  = 2, jmax1
  3903.           ioff = (j-1) * imaxl + ibeg
  3904.        do 10 i = 2, imax1
  3905.             ij = i + ioff
  3906.         sum    = yfu(ij) + yco2(ij) + yh2o(ij) + yo2(ij) + yn2(ij)
  3907.         yfu (ij) = yfu (ij) / sum
  3908.         yco2(ij) = yco2(ij) / sum
  3909.         yh2o(ij) = yh2o(ij) / sum
  3910.         yo2 (ij) = yo2 (ij) / sum
  3911.         yn2 (ij) = yn2 (ij) / sum
  3912. 10     continue
  3913. 11    continue
  3914. C
  3915.       return
  3916.       end
  3917. C======================================================================C
  3918.       subroutine trsrc (igrl, q)
  3919. C                                                                      C
  3920. C     Purpose:  Compute the source terms resulting from coordinate     C
  3921. C               transformation.                                        C
  3922. C                                                                      C
  3923. C======================================================================C
  3924.       include 'UIFlow.com'
  3925.       dimension q(*)
  3926.       include 'UIFlow.indx'
  3927. C
  3928. C.... Compute the derivatives at the cell centers.
  3929. C
  3930.       do 11 j  = 2, jmax1
  3931.           ioff = (j-1) * imaxl + ibeg
  3932.          ioff1 = ioff - 1
  3933.          ioffm = ioff - imaxl
  3934.        do 10 i = 2, imax1
  3935.            ij  = i + ioff
  3936.           ij1  = i + ioff1
  3937.           ijm  = i + ioffm
  3938.         x1(ij) = fx(ij)  * q(ij+1) + fx1(ij)  * q(ij) -
  3939.      1           fx(ij1) * q(ij)   - fx1(ij1) * q(ij1)
  3940.         y1(ij) = fy(ij)  * q(ij+imaxl) + fy1(ij) * q(ij) -
  3941.      1           fy(ijm) * q(ij)   - fy1(ijm) * q(ijm)
  3942.         x1(ij) = gam(ij) * q21(ij) * x1(ij)
  3943.         y1(ij) = gam(ij) * q12(ij) * y1(ij)
  3944. 10     continue
  3945. 11    continue
  3946. C
  3947. C.... Calculate derivative at the x - minus boundary.
  3948. C
  3949.       do 20 j = 2, jmax1
  3950.          ij   = 1 + (j-1) * imaxl + ibeg
  3951.         ijm   = ij - imaxl
  3952.        y1(ij) = fy(ij)  * q(ij+imaxl) + fy1(ij)  * q(ij) -
  3953.      1          fy(ijm) * q(ij)       - fy1(ijm) * q(ijm)
  3954.        y1(ij) = gam(ij) * q12(ij+1) * y1(ij)
  3955. 20    continue
  3956. C
  3957. C.... Calculate derivative at the y - minus boundary.
  3958. C
  3959.       do 21 i = 2, imax1
  3960.          ij   = ibeg + i
  3961.          ij1  = ij + imaxl
  3962.          ijm  = ij - 1
  3963.        x1(ij) = fx(ij)  * q(ij+1) + fx1(ij)  * q(ij) -
  3964.      1          fx(ijm) * q(ij)   - fx1(ijm) * q(ijm)
  3965.        x1(ij) = gam(ij) * q21(ij1) * x1(ij)
  3966. 21    continue
  3967. C
  3968. C.... Calculate derivative at the x - plus boundary.
  3969. C
  3970.       do 22 j = 2, jmax1
  3971.           ij  = imaxl + (j-1) * imaxl + ibeg
  3972.          ij1  = ij - 1
  3973.          ijm  = ij - imaxl
  3974.        y1(ij) = fy(ij)  * q(ij+imaxl) + fy1(ij)  * q(ij) -
  3975.      1          fy(ijm) * q(ij)       - fy1(ijm) * q(ijm)
  3976.        y1(ij) = gam(ij) * q12(ij1) * y1(ij)
  3977. 22    continue
  3978. C
  3979. C.... Calculate derivative at the y - plus boundary.
  3980. C
  3981.       do 23 i = 2, imax1
  3982.          ij   = i + jmax1 * imaxl + ibeg 
  3983.          ij1  = ij - imaxl
  3984.          ijm  = ij - 1
  3985.        x1(ij) = fx(ij)  * q(ij+1) + fx1(ij)  * q(ij) -
  3986.      1          fx(ijm) * q(ij)   - fx1(ijm) * q(ijm)
  3987.        x1(ij) = gam(ij)  * q21(ij1) * x1(ij)
  3988. 23    continue
  3989. C
  3990. C.... Calculate transformation source term.
  3991. C
  3992.       do 51 j   = 2, jmax1
  3993.            ioff = (j-1) * imaxl + ibeg
  3994.        do 50 i  = 2, imax1
  3995.             ij  = i + ioff
  3996.             ij1 = ij - 1
  3997.             ijp = ij + imaxl
  3998.             ijm = ij - imaxl
  3999.         su(ij)  = ( fx(ij)  * y1(ij+1) + fx1(ij)  * y1(ij)   )
  4000.      1          - ( fx(ij1) * y1(ij)   + fx1(ij1) * y1(ij1) )
  4001.      2          + ( fy(ij)  * x1(ijp)  + fy1(ij)  * x1(ij)   )
  4002.      3          - ( fy(ijm) * x1(ij)   + fy1(ijm) * x1(ijm)  )
  4003. 50     continue
  4004. 51    continue
  4005.       return
  4006.       end
  4007. C======================================================================C
  4008.       subroutine tvis (igrl)
  4009. C                                                                      C
  4010. C     Purpose:  Calculate the turbulent viscosity.  Together with      C
  4011. C               the laminar viscosity it comprises the exchange        C
  4012. C               coefficient for the momentum equations.                C
  4013. C                                                                      C
  4014. C======================================================================C
  4015.       include 'UIFlow.com'
  4016.       include 'UIFlow.indx'
  4017. C
  4018.       relxm    = 1.0 - relx(6)
  4019. C
  4020.       do 11 j  = 2, jmax1
  4021.          ioff  = ibeg + (j-1)* imaxl
  4022.        do 10 i = 2, imax1
  4023.            ij  = i + ioff
  4024.         amut(ij) = relxm * amut(ij) + 
  4025.      1           relx(6) * cd * rho(ij) * tke(ij) * tke(ij) / tde(ij)
  4026. 10     continue
  4027. 11    continue
  4028. C
  4029.       return
  4030.       end
  4031. C======================================================================C
  4032.       subroutine update (igrf, igrl)
  4033. C                                                                      C
  4034. C     Purpose:  Correct the flow variables in accordance with the use  C
  4035. C               of a decoupled solution procedure ( SIMPLEC ).         C
  4036. C                                                                      C
  4037. C======================================================================C
  4038.       include 'UIFlow.com'
  4039.       include 'UIFlow.indx'
  4040. C
  4041.       call bcor (igrl, pp)
  4042.       call grad (igrl, pp)
  4043. C
  4044. C.... Corrections for incompressible flows.
  4045. C
  4046.       do 11 j  = 2, jmax1
  4047.          ioff  = (j-1) * imaxl + ibeg
  4048.        do 10 i = 2, imax1
  4049.            ij  = i + ioff
  4050.         rapm   = 1.0 / apm(ij)
  4051.         p(ij)  = p(ij) +  pp(ij) * relx(3)
  4052.         u(ij)  = u(ij) + (dux(ij) * dpdx(ij) + duy(ij) * dpdy(ij))
  4053.      1                 * rapm
  4054.         v(ij)  = v(ij) + (dvx(ij) * dpdx(ij) + dvy(ij) * dpdy(ij))
  4055.      1                 * rapm
  4056. 10     continue
  4057. 11    continue
  4058. C
  4059. C.... Correct zi - direction fluxes.
  4060. C
  4061.       do 13 j  = 2, jmax1
  4062.          ioff  = (j-1) * imaxl + ibeg
  4063.        do 12 i = 2, imax2
  4064.            ij  = i + ioff
  4065.           ij1  = ij + 1
  4066.         cx(ij) = cx(ij) + dx(ij) * ( pp(ij) - pp(ij1) )
  4067.         cu(ij) = cx(ij) / ( amax1( sign( rho(ij) ,  cx(ij) ), 0.)
  4068.      1                  +   amax1( sign( rho(ij1), -cx(ij) ), 0.) )
  4069. 12     continue
  4070. 13    continue
  4071. C
  4072. C.... Correct eta direction fluxes.
  4073. C
  4074.       do 15 j = 2, jmax2
  4075.          ioff = (j-1) * imaxl + ibeg
  4076.         ioffp = ioff + imaxl
  4077.        do 14 i = 2, imax1
  4078.            ij = i + ioff
  4079.           ijp = i + ioffp
  4080.         cy(ij)  = cy(ij) + dy(ij) * ( pp(ij) - pp(ijp) )
  4081.         cv(ij)  = cy(ij) / ( amax1( sign( rho(ij) ,  cv(ij) ), 0.)
  4082.      1                   +   amax1( sign( rho(ijp), -cv(ij) ), 0.) )
  4083. 14     continue
  4084. 15    continue
  4085. C
  4086. C.... Additional corrections for compressibility.
  4087. C
  4088.       if ( kcomp .eq. 1 ) then
  4089.        do 21 j  = 2, jmax1
  4090.           ioff  = (j-1) * imaxl + ibeg
  4091.          ioffp  = ioff + imaxl
  4092.         do 20 i = 2, imax2
  4093.             ij  = i + ioff
  4094.            ijp  = i + ioffp
  4095.            ij1  = ij + 1
  4096.          cx(ij) = cx(ij) + amax1( -cu(ij), 0.) * wmol(ij1) * pp(ij1)
  4097.      1                         / (  gamma * gascon * t(ij1) )
  4098.      2                   + amax1(  cu(ij), 0.) * wmol(ij) * pp(ij)
  4099.      3                         / (  gamma * gascon * t(ij)   )
  4100. 20      continue
  4101. 21     continue
  4102. C
  4103.        do 23 j  = 2, jmax2
  4104.           ioff  = (j-1) * imaxl + ibeg
  4105.          ioffp  = ioff + imaxl
  4106.         do 22 i = 2, imax1
  4107.             ij  = i + ioff
  4108.            ijp  = i + ioffp
  4109.          cy(ij) = cy(ij) + amax1( -cv(ij), 0. ) * wmol(ijp) * pp(ijp)
  4110.      1                       / (  gamma * gascon * t(ijp) )
  4111.      2                   + amax1(  cv(ij), 0.) * wmol(ij) * pp(ij)
  4112.      3                       / (  gamma * gascon * t(ij)  )
  4113. 22      continue
  4114. 23     continue
  4115.       endif
  4116. C
  4117. C.... Corrections to fluxes for highly non-orthogonal grid.
  4118. C
  4119.       if ( knorth .eq. 1 ) then
  4120.        do 32 j  = 2, jmax1
  4121.           ioff  = (j-1) * imaxl + ibeg
  4122.         do 31 i = 2, imax2
  4123.              ij = i + ioff
  4124.             ij1 = ij + 1
  4125.             ijp = ij + imaxl
  4126.          cu(ij) = a11(ij) * ( fx1(ij)*duy(ij) *dpdy(ij) /apm(ij)
  4127.      1                    +   fx(ij) *duy(ij1)*dpdy(ij1)/apm(ij1) )
  4128.      2          + a12(ij) * ( fx1(ij)*dvy(ij) *dpdy(ij) /apm(ij)
  4129.      3                    +   fx(ij) *dvy(ij1)*dpdy(ij1)/apm(ij1) )
  4130.      4          + cu(ij)
  4131.          cx(ij) = cu(ij) * ( amax1( sign( rho(ij) ,  cu(ij) ), 0.)
  4132.      1                   +   amax1( sign( rho(ij1), -cu(ij) ), 0.) )
  4133. 31      continue
  4134. 32     continue
  4135. C
  4136.        do 34 j  = 2, jmax2
  4137.           ioff  = (j-1) * imaxl + ibeg
  4138.         do 33 i = 2, imax1
  4139.             ij  = i + ioff
  4140.            ijp  = ij + imaxl
  4141.          cv(ij) = a21(ij) * ( fy1(ij)*dux(ij) *dpdx(ij) /apm(ij)
  4142.      1                    +   fy(ij) *dux(ijp)*dpdx(ijp)/apm(ijp) )
  4143.      2          + a22(ij) * ( fy1(ij)*dvx(ij) *dpdx(ij) /apm(ij)
  4144.      3                    +   fy(ij) *dvx(ijp)*dpdx(ijp)/apm(ijp) )
  4145.      4          + cv(ij)
  4146.          cy(ij) = cv(ij)  * ( amax1( sign( rho(ij) ,  cv(ij) ), 0.)
  4147.      1                    +   amax1( sign( rho(ijp), -cv(ij) ), 0.) )
  4148. 33      continue
  4149. 34     continue
  4150.       endif
  4151. C
  4152.       return
  4153.       end
  4154. C======================================================================C
  4155.       subroutine urelax(igrl,nv,q)
  4156. C                                                                      C
  4157. C     Purpose:  Incorporate under-relaxation into the discretized      C
  4158. C               equation.                                              C
  4159. C                                                                      C
  4160. C======================================================================C
  4161.       include 'UIFlow.com'
  4162.       dimension q(*)
  4163.       include 'UIFlow.indx'
  4164. C
  4165.       relxm   = 1.0 - relx(nv)
  4166. C
  4167.       do 11 j  = 2, jmax1
  4168.          ioff  = (j-1) * imaxl + ibeg
  4169.        do 10 i = 2, imax1
  4170.            ij  = i + ioff
  4171.         ap(ij) = ae(ij) + aw(ij) + an(ij) + as(ij) - sp(ij)
  4172.         sp(ij) = 0.0
  4173.         ap(ij) = ap(ij) / relx(nv)
  4174.         su(ij) = su(ij) + relxm * ap(ij) * q(ij)
  4175. 10     continue
  4176. 11    continue
  4177. C
  4178.       if ( nv .eq. 1 ) then
  4179.        do 20  j  = 2, jmax1
  4180.            ioff  = (j-1) * imaxl + ibeg
  4181.         do 21 i  = 2, imax1
  4182.              ij  = i + ioff
  4183.          app(ij) = ap(ij)
  4184. 21      continue
  4185. 20     continue
  4186.       endif
  4187. C
  4188.       return
  4189.       end
  4190. C======================================================================C
  4191.       subroutine visc (igrl, nv)
  4192. C                                                                      C
  4193. C     Purpose:  Calculate diffusive exchange coefficient.              C
  4194. C                                                                      C
  4195. C======================================================================C
  4196.       include 'UIFlow.com'
  4197.       include 'UIFlow.indx'
  4198. C
  4199.       do 11 j   = 1, jmaxl
  4200.           ioff  = (j-1) * imaxl + ibeg
  4201.        do 10 i  = 1, imaxl
  4202.             ij  = i + ioff
  4203.         gam(ij) = amu(ij) / prl(nv) + amut(ij) / prt(nv)
  4204. 10     continue
  4205. 11    continue
  4206. C
  4207.       return
  4208.       end
  4209. C======================================================================C
  4210.       subroutine vset
  4211. C======================================================================C
  4212.       include 'UIFlow.com'
  4213. C
  4214.       dimension xloc(3500) , yloc(3500) , uf(3500)   , vf(3500)  ,
  4215.      1          wf  (3500) , tf(3500)   , ff(3500)   , gf(3500)  ,
  4216.      2          yfuf(3500) , yh2of(3500), yco2f(3500), yo2f(3500),
  4217.      3          yn2f(3500) , rhof(3500) , hf(3500)   , tkef(3500),
  4218.      4          tdef(3500) , pf(3500)   , n1(14000)
  4219.       external DFopen, DFclose, vsfatch, vsfdtch, vsfsfld
  4220.       external vsfsnam, vsfwrit, vsffdef, vfatch, vfdtch, vfsnam
  4221.       external vfinsrt
  4222.       integer  DFopen, vsfatch, vsfsfld, vsfwrit, vsffdef, vfatch,
  4223.      1         vfinsrt, ftype, fintrlace, fullacc,
  4224.      2         vg, vs
  4225. C    ==================================================================
  4226. C
  4227.       external VHFSD, VHFSDM, VHFMKGP
  4228.       integer    VHFSD, VHFSDM, VHFMKGP
  4229.       integer INTTYPE    
  4230.       character*10 strng
  4231.       character*20 strng1
  4232.       character*20 strng2
  4233.       character*10 uiflow
  4234.       parameter (INTTYPE=2)
  4235.       integer REALTYPE                
  4236.       parameter (REALTYPE=3)
  4237.       integer    VDATATAG
  4238.       parameter (VDATATAG=1962)
  4239.       integer tagarray(30), refarray(30), vxsize
  4240. C    ==================================================================
  4241. C    
  4242.       CEXTERNAL SHOWLINE
  4243.       CHARACTER*255 buffer
  4244.       CHARACTER*1 NULL
  4245. C
  4246.       parameter (ftype = 3, fintrlace = 0, fullacc = 7 )
  4247.       igrl = ngrid
  4248.       include 'UIFlow.indx'
  4249.       NULL = char(0)
  4250. C
  4251.       nsize = imaxl * jmaxl
  4252.       nsize2 = imax1 * jmax1
  4253. C
  4254. C.... Calculate x-y coordinates at data points.
  4255. C
  4256.       do 10 j  = 2, jmax1
  4257.        do 11 i = 2, imax1
  4258.             ij = ibeg + i + (j-1) * imaxl
  4259.            ijm = ij - imaxl
  4260.           i0j0 = i + (j-1) * imaxl
  4261.         xloc(i0j0) = 0.25 * ( x(ij) + x(ij-1) + x(ijm) + x(ijm-1) )
  4262.         yloc(i0j0) = 0.25 * ( y(ij) + y(ij-1) + y(ijm) + y(ijm-1) )
  4263. 11     continue
  4264. 10    continue
  4265. C
  4266. C.... Calculate y-minus and y-plus x,y locations.
  4267. C
  4268.       do 12 i = 2, imax1
  4269.           ijm = ibeg + i
  4270.           ijp = ibeg + i + (jmax2) * imaxl
  4271.          i0jm = i
  4272.          i0jp = i + (jmax1) * imaxl
  4273.        xloc(i0jm) = 0.5 * ( x(ijm) + x(ijm-1) )
  4274.        xloc(i0jp) = 0.5 * ( x(ijp) + x(ijp-1) )
  4275.        yloc(i0jm) = 0.5 * ( y(ijm) + y(ijm-1) )
  4276.        yloc(i0jp) = 0.5 * ( y(ijp) + y(ijp-1) )
  4277. 12    continue
  4278. C
  4279. C.... Calculate x-minus and x-plus x,y locations.
  4280. C
  4281.       do 13 j = 2, jmax1
  4282.           imj = ibeg + 1 + (j-1) * imaxl
  4283.           ipj = ibeg + imax1 + (j-1) * imaxl
  4284.          imj0 = 1 + (j-1) * imaxl
  4285.          ipj0 = imaxl + (j-1) * imaxl
  4286.        xloc(imj0) = 0.5 * ( x(imj) + x(imj-imaxl) )
  4287.        xloc(ipj0) = 0.5 * ( x(ipj) + x(ipj-imaxl) )
  4288.        yloc(imj0) = 0.5 * ( y(imj) + y(imj-imaxl) )
  4289.        yloc(ipj0) = 0.5 * ( y(ipj) + y(ipj-imaxl) )
  4290. 13    continue
  4291. C
  4292. C.... Specify CORNER points.
  4293. C
  4294.       xloc(1) = x(ibeg + 1)
  4295.       yloc(1) = y(ibeg + 1)
  4296. C
  4297.       xloc(imaxl) = x(ibeg + imax1)
  4298.       yloc(imaxl) = y(ibeg + imax1)
  4299. C
  4300.       xloc(1 + jmax1*imaxl) = x(ibeg + 1 + jmax2*imaxl)
  4301.       yloc(1 + jmax1*imaxl) = y(ibeg + 1 + jmax2*imaxl)
  4302. C
  4303.       xloc(imaxl + jmax1*imaxl) = x(ibeg + imax1 + jmax2*imaxl)
  4304.       yloc(imaxl + jmax1*imaxl) = y(ibeg + imax1 + jmax2*imaxl)
  4305. C
  4306.       do 20 j   = 1, jmaxl
  4307.        do 21 i  = 1, imaxl
  4308.             ij  = ibeg + i + (j-1) * imaxl
  4309.           i0j0  = i + (j-1 ) * imaxl
  4310.         uf   (i0j0) = u   (ij)
  4311.         vf   (i0j0) = v   (ij)
  4312.         wf   (i0j0) = w   (ij)
  4313.         pf   (i0j0) = p   (ij)
  4314.         tf   (i0j0) = t   (ij)
  4315.         hf   (i0j0) = h   (ij)
  4316.         ff   (i0j0) = f   (ij)
  4317.         gf   (i0j0) = g   (ij)
  4318.         rhof (i0j0) = rho (ij)
  4319.         yfuf (i0j0) = yfu (ij)
  4320.         yh2of(i0j0) = yh2o(ij)
  4321.         yco2f(i0j0) = yco2(ij)
  4322.         yo2f (i0j0) = yo2 (ij)
  4323.         yn2f (i0j0) = yn2 (ij)
  4324.         tkef (i0j0) = tke (ij)
  4325.         tdef (i0j0) = tde (ij)
  4326. 21     continue
  4327. 20    continue
  4328. C
  4329. C.... Compute connectivity.
  4330. C
  4331.       do 30 j  = 1, jmax1
  4332.        do 31 i = 1, imax1
  4333.             ij = i + (j-1) * imaxl
  4334.            ncv = i + (j-1) * imax1
  4335.           indx = ( ncv - 1 ) * 4 + 1
  4336.         n1(indx)   = ij
  4337.         n1(indx+1) = ij + 1
  4338.         n1(indx+3) = ij + imaxl
  4339.         n1(indx+2) = ij + imaxl + 1
  4340. 31     continue
  4341. 30    continue
  4342. C
  4343. C
  4344.       nn = DFopen('vset.out', fullacc, 0 )
  4345.       if (nn .eq. -1) then 
  4346.           call SHOWLINE('Unable to Create VSet')
  4347.         return
  4348.       endif
  4349. C
  4350.       vxsize = 0
  4351. C
  4352.       uiflow = 'uiflow' // char(0)
  4353.       strng = 'px' // char(0)
  4354.       strng1 = 'xcoord' // char(0)
  4355.       vs = VHFSD (nn, strng, xloc, nsize, REALTYPE, strng1, uiflow)
  4356.       vxsize = vxsize + 1
  4357.       tagarray(vxsize) = VDATATAG
  4358.       refarray(vxsize) = vs
  4359. C
  4360.       strng = 'py' // char(0)
  4361.       strng1 = 'ycoord' // char(0)
  4362.       vs = VHFSD (nn, strng, yloc, nsize, REALTYPE, strng1, uiflow)
  4363.       vxsize = vxsize + 1
  4364.       tagarray(vxsize) = VDATATAG
  4365.       refarray(vxsize) = vs
  4366. C
  4367.       strng = 'plist4' // char(0)
  4368.       strng1 = 'CONNECTIVITY' // char(0)
  4369.       vs = VHFSDM (nn, strng, n1, nsize2, INTTYPE, strng1,uiflow,4)
  4370.       vxsize = vxsize + 1
  4371.       tagarray(vxsize) = VDATATAG
  4372.       refarray(vxsize) = vs
  4373. C
  4374. C.... Write data.
  4375. C
  4376.       if (iprint(1) .eq. 1) then
  4377.         strng = 'u-vel' // char(0)
  4378.         strng1 = 'U-VELOCITY' // char(0)
  4379.         vs = VHFSD (nn, strng, uf, nsize, REALTYPE, strng1,uiflow)
  4380.         vxsize = vxsize + 1
  4381.         tagarray(vxsize) = VDATATAG
  4382.         refarray(vxsize) = vs
  4383.       endif
  4384. C
  4385.       if (iprint(2) .eq. 1) then
  4386.         strng = 'v-vel' // char(0)
  4387.         strng1 = 'V-VELOCITY' // char(0)
  4388.         vs = VHFSD (nn, strng, vf, nsize, REALTYPE, strng1,uiflow)
  4389.         vxsize = vxsize + 1
  4390.         tagarray(vxsize) = VDATATAG
  4391.         refarray(vxsize) = vs
  4392.       endif
  4393. C
  4394.       if (iprint(4) .eq. 1) then
  4395.         strng = 'w-vel' // char(0)
  4396.         strng1 = 'W-VELOCITY' // char(0)
  4397.         vs = VHFSD (nn, strng, wf, nsize, REALTYPE, strng1,uiflow)
  4398.         vxsize = vxsize + 1
  4399.         tagarray(vxsize) = VDATATAG
  4400.         refarray(vxsize) = vs
  4401.       endif
  4402. C
  4403.       if (iprint(5) .eq. 1) then
  4404.         strng = 'temp' // char(0)
  4405.         strng1 = 'TEMPERATURE' // char(0)  
  4406.         vs = VHFSD (nn, strng, tf, nsize, REALTYPE, strng1,uiflow)
  4407.         vxsize = vxsize + 1
  4408.         tagarray(vxsize) = VDATATAG
  4409.         refarray(vxsize) = vs
  4410.       endif
  4411. C
  4412.       if (iprint(3) .eq. 1) then
  4413.         strng = 'press' // char(0)
  4414.         strng1 = 'PRESSURE' // char(0)
  4415.         vs = VHFSD (nn, strng, pf, nsize, REALTYPE, strng1,uiflow)
  4416.         vxsize = vxsize + 1
  4417.         tagarray(vxsize) = VDATATAG
  4418.         refarray(vxsize) = vs
  4419.       endif
  4420. C
  4421.       if (iprint(11) .eq. 1) then
  4422.         strng = 'dens' // char(0)
  4423.         strng1 = 'DENSITY' // char(0)
  4424.         vs = VHFSD (nn, strng, rhof, nsize, REALTYPE, strng1,uiflow)
  4425.         vxsize = vxsize + 1
  4426.         tagarray(vxsize) = VDATATAG
  4427.         refarray(vxsize) = vs
  4428.       endif
  4429. C
  4430.       if (iprint(5) .eq. 1) then
  4431.         strng = 'enth' // char(0)
  4432.         strng1 = 'ENTHALPY' // char(0)
  4433.         vs = VHFSD (nn, strng, hf, nsize, REALTYPE, strng1,uiflow)
  4434.         vxsize = vxsize + 1
  4435.         tagarray(vxsize) = VDATATAG
  4436.         refarray(vxsize) = vs
  4437.       endif
  4438. C
  4439.       if (iprint(9) .eq. 1) then
  4440.         strng = 'fuel' // char(0)
  4441.         strng1 = 'FUEL' // char(0)
  4442.         vs = VHFSD (nn,strng, yfuf, nsize, REALTYPE, strng1,uiflow)
  4443.         vxsize = vxsize + 1
  4444.         tagarray(vxsize) = VDATATAG
  4445.         refarray(vxsize) = vs
  4446.       endif
  4447. C
  4448.       if (iprint(12) .eq. 1) then
  4449.         strng = 'h2o' // char(0)
  4450.         strng1 = 'WATER' // char(0)
  4451.         vs = VHFSD (nn, strng, yh2of, nsize, REALTYPE, strng1,uiflow)
  4452.         vxsize = vxsize + 1
  4453.         tagarray(vxsize) = VDATATAG
  4454.         refarray(vxsize) = vs
  4455.       endif
  4456. C
  4457.       if (iprint(12) .eq. 1) then
  4458.         strng = 'co2' // char(0)
  4459.         strng1 = 'CARBON_DIOXIDE' // char(0)
  4460.         vs = VHFSD (nn, strng, yco2f, nsize, REALTYPE, strng1,uiflow)
  4461.         vxsize = vxsize + 1
  4462.         tagarray(vxsize) = VDATATAG
  4463.         refarray(vxsize) = vs
  4464.       endif
  4465. C
  4466.       if (iprint(12) .eq. 1) then
  4467.         strng = 'o2' // char(0)
  4468.         strng1 = 'OXYGEN' // char(0)
  4469.         vs = VHFSD (nn, strng, yo2f, nsize, REALTYPE, strng1,uiflow)
  4470.         vxsize = vxsize + 1
  4471.         tagarray(vxsize) = VDATATAG
  4472.         refarray(vxsize) = vs
  4473.       endif
  4474. C
  4475.       if (iprint(12) .eq. 1) then
  4476.         strng = 'n2' // char(0)
  4477.         strng1 = 'NITROGEN' // char(0)
  4478.         vs = VHFSD (nn, strng, yn2f, nsize, REALTYPE, strng1,uiflow)
  4479.         vxsize = vxsize + 1
  4480.         tagarray(vxsize) = VDATATAG
  4481.         refarray(vxsize) = vs
  4482.       endif
  4483. C
  4484.       if (iprint(8) .eq. 1) then
  4485.         strng = 'mix.frac.' // char(0)
  4486.         strng1 = 'MIXTURE.FRAC' // char(0)
  4487.         vs = VHFSD (nn, strng, ff, nsize, REALTYPE, strng1,uiflow)
  4488.         vxsize = vxsize + 1
  4489.         tagarray(vxsize) = VDATATAG
  4490.         refarray(vxsize) = vs
  4491.       endif
  4492. C
  4493.       if (iprint(10) .eq. 1) then
  4494.         strng = 'conc.fluc.' // char(0)
  4495.         strng1 = 'CONC.FLUC.' // char(0)
  4496.         vs = VHFSD (nn, strng, gf, nsize, REALTYPE, strng1, uiflow)
  4497.         vxsize = vxsize + 1
  4498.         tagarray(vxsize) = VDATATAG
  4499.         refarray(vxsize) = vs
  4500.       endif
  4501. C
  4502.       if (iprint(6) .eq. 1) then
  4503.         strng = 'kin.ener.' // char(0)
  4504.         strng1 = 'KIN.ENER.' // char(0)
  4505.         vs = VHFSD (nn, strng, tkef, nsize, REALTYPE, strng1,uiflow)
  4506.         vxsize = vxsize + 1
  4507.         tagarray(vxsize) = VDATATAG
  4508.         refarray(vxsize) = vs
  4509.       endif
  4510. C
  4511.       if (iprint(7) .eq. 1) then
  4512.         strng = 'turb.disp' // char(0)
  4513.         strng1 = 'TURB.DISP.' // char(0)
  4514.         vs = VHFSD (nn, strng, tdef, nsize, REALTYPE, strng1,uiflow)
  4515.         vxsize = vxsize + 1
  4516.         tagarray(vxsize) = VDATATAG
  4517.         refarray(vxsize) = vs
  4518.       endif
  4519. C
  4520.       strng2 = 'UIFlow Group' // char(0)
  4521.       vg = VHFMKGP(nn,tagarray,refarray,vxsize,strng2,uiflow)
  4522.       if (vg .eq. -1)   call SHOWLINE('Unable to Create VGroup')
  4523.       call DFclose ( nn )
  4524. C
  4525.       return
  4526.       end
  4527. C======================================================================C
  4528.       subroutine walld (igrl,i1,i2,j1,j2,dst)
  4529. C                                                                      C
  4530. C     Purpose:  Modify the source terms at the walls to impose wall    C
  4531. C               functions on turbulent dissipation.                    C
  4532. C                                                                      C
  4533. C======================================================================C
  4534.       include 'UIFlow.com'
  4535.       dimension dst(*)
  4536.       include 'UIFlow.indx'
  4537. C
  4538. C.... Turb. dissipation is fixed at near wall nodes.
  4539. C
  4540.       do 11 i  = i1, i2
  4541.        do 10 j = j1, j2
  4542.            ij  = i + (j-1) * imaxl + ibeg
  4543.         su(ij) = 2.0e10 * cdtqtr * (tke(ij) ** 1.5) / (ak * dst(ij))
  4544.         sp(ij) = - 1.0e10
  4545. 10     continue
  4546. 11    continue
  4547. C
  4548.       return
  4549.       end
  4550. C======================================================================C
  4551.       subroutine wallg (igrl, i1, i2, j1, j2, ioff, dst)
  4552. C                                                                      C
  4553. C     Purpose:  Modify the values of the exchange coefficient at the   C
  4554. C               walls so as to impose wall functions on velocities.    C
  4555. C                                                                      C
  4556. C======================================================================C
  4557.       include 'UIFlow.com'
  4558.       dimension dst(*)
  4559.       include 'UIFlow.indx'
  4560. C
  4561.       do 11 i  = i1, i2
  4562.        do 10 j = j1, j2
  4563.            ij  = i + (j-1) * imaxl + ibeg
  4564.           ij1  = ij + ioff
  4565.         term     = 0.5  * rho(ij) * cdqtr * sqrt( tke(ij) ) * dst(ij)
  4566.         gam(ij1) = term * ak / alog (ee * term / amu(ij) )
  4567. 10     continue
  4568. 11    continue
  4569. C
  4570.       return
  4571.       end
  4572. C======================================================================C
  4573.       subroutine wallk(igrl,i1,i2,j1,j2,v1,v2,v3,vw1,vw2,vw3,dst,cf1)
  4574. C                                                                      C
  4575. C     Purpose:  Modify the source terms at the near wall cells so as   C
  4576. C               to impose wall functions on k.                         C
  4577. C                                                                      C
  4578. C======================================================================C
  4579.       include 'UIFlow.com'
  4580.       dimension v1(*), v2(*), v3(*), cf1(*), dst(*)
  4581.       include 'UIFlow.indx'
  4582. C
  4583.       do 11 j   = j1, j2
  4584.        do 10 i  = i1, i2
  4585.             ij  = i + (j-1) * imaxl + ibeg
  4586.         tkhf    =   sqrt ( tke(ij) )
  4587.         vres    =   sqrt ( u(ij)*u(ij) + v(ij)*v(ij) + w(ij)*w(ij) )
  4588.         dist    =   0.5 * dst(ij)
  4589.         term    =   rho(ij) * cdqtr * tkhf
  4590.         denom   =   alog ( ee * dist * term / amu(ij) )
  4591.         tau     =   ak * term * vres / denom
  4592.         tau1    =   amu(ij) * vres  / dist
  4593.         if (tau .lt. 0.0) tau = tau1
  4594.         prdn    =   tau * vres * ajb(ij) / dist
  4595.         sp(ij)  = - rho(ij) * cdtqtr * tkhf * denom / 
  4596.      1              (ak * dist) * ajb(ij)
  4597.         su(ij)  =   prdn
  4598.         cf1(ij) =   0.0
  4599. 10     continue
  4600. 11    continue
  4601. C
  4602.       return
  4603.       end
  4604. C======================================================================C
  4605.       subroutine xmextr (igrl)
  4606. C                                                                      C
  4607. C     Purpose:  Extrapolate appropriate boundary conditions to the     C
  4608. C               x-minus ( zi-minus ) boundary.                         C
  4609. C                                                                      C
  4610. C======================================================================C
  4611.       include 'UIFlow.com'
  4612.       include 'UIFlow.indx'
  4613. C
  4614.       i = 1
  4615. C
  4616.       do 100 n = 1, nsxm
  4617.           jfl  = jfxm(n,igrl)
  4618.       jll  = jlxm(n,igrl)
  4619. C
  4620. C.... Enforce WALL boundary conditions.
  4621. C
  4622.       if ( kbxm(n) .eq. 1 ) then
  4623.        do 10 j   = jfl, jll
  4624.             ij   = ibeg + i + (j-1) * imaxl
  4625.            ij1   = ij + 1
  4626.         p   (ij) = p   (ij1)
  4627.         rho (ij) = rho (ij1)
  4628.         tke (ij) = tke (ij1)
  4629.         tde (ij) = tde (ij1)
  4630.         amut(ij) = amut(ij1)
  4631.         h   (ij) = h   (ij1)
  4632.         t   (ij) = t   (ij1)
  4633.         f   (ij) = f   (ij1)
  4634.         g   (ij) = g   (ij1)
  4635.         yfu (ij) = yfu (ij1)
  4636.         yo2 (ij) = yo2 (ij1)
  4637.         yco2(ij) = yco2(ij1)
  4638.         yh2o(ij) = yh2o(ij1)
  4639.         yn2 (ij) = yn2 (ij1)
  4640. 10     continue
  4641. C
  4642. C.... Enforce INLET boundary condition.
  4643. C
  4644.       elseif ( kbxm(n) .eq. 2 ) then
  4645.        do 20 j = jfl, jll
  4646.            ij  = ibeg + i + (j-1) * imaxl
  4647.            ij1 = ij + 1
  4648.         p(ij)  = p (ij1) + 0.5 * ( p(ij1) - p(ij+2) )
  4649. 20     continue
  4650. C
  4651. C.... Enforce SYMMETRY boundary condition.
  4652. C
  4653.       elseif ( kbxm(n) .eq. 3 ) then
  4654.        do 30 j = jfl, jll
  4655.            ij  = ibeg + i + (j-1) * imaxl
  4656.            ij1 = ij + 1
  4657.         cx   (ij) = 0.0
  4658.         cu   (ij) = 0.0
  4659.         cy   (ij) = cy (ij1)
  4660.         rho  (ij) = rho(ij1)
  4661.         cv   (ij) = cy (ij) / rho(ij)
  4662.         p    (ij) = p   (ij1)
  4663.         tke  (ij) = tke (ij1)
  4664.         tde  (ij) = tde (ij1)
  4665.         amut (ij) = amut(ij1)
  4666.         h    (ij) = h   (ij1)
  4667.         t    (ij) = t   (ij1)
  4668.         u    (ij) = u   (ij1)
  4669.         v    (ij) = v   (ij1)
  4670.         w    (ij) = w   (ij1)
  4671.         f    (ij) = f   (ij1)
  4672.         g    (ij) = g   (ij1)
  4673.         yfu  (ij) = yfu (ij1)
  4674.         yo2  (ij) = yo2 (ij1)
  4675.         yco2 (ij) = yco2(ij1)
  4676.         yh2o (ij) = yh2o(ij1)
  4677.         yn2  (ij) = yn2 (ij1)
  4678. 30     continue
  4679. C
  4680. C.... Enforce OUTFLOW boundary condition.
  4681. C
  4682.       elseif ( kbxm(n) .eq. 4 ) then
  4683.        do 40 j = jfl, jll
  4684.            ij  = ibeg + i + (j-1) * imaxl
  4685.            ij1 = ij + 1
  4686.         cx   (ij) = cx (ij1)
  4687.         rho  (ij) = rho(ij1)
  4688.         cu   (ij) = cx(ij) / rho(ij)
  4689.         cy   (ij) = 0.0
  4690.         cv   (ij) = 0.0
  4691.         p    (ij) = p (ij1) + 0.5 * ( p(ij1) - p(ij+2) )
  4692.         tke  (ij) = tke (ij1)
  4693.         tde  (ij) = tde (ij1)
  4694.         amut (ij) = amut(ij1)
  4695.         h    (ij) = h   (ij1)
  4696.         t    (ij) = t   (ij1)
  4697.         u    (ij) = u   (ij1)
  4698.         v    (ij) = v   (ij1)
  4699.         w    (ij) = w   (ij1)
  4700.         f    (ij) = f   (ij1)
  4701.         g    (ij) = g   (ij1)
  4702.         yfu  (ij) = yfu (ij1)
  4703.         yo2  (ij) = yo2 (ij1)
  4704.         yco2 (ij) = yco2(ij1)
  4705.         yh2o (ij) = yh2o(ij1)
  4706.         yn2  (ij) = yn2 (ij1)
  4707. 40     continue
  4708.       endif
  4709. C
  4710. 100   continue
  4711.       return
  4712.       end
  4713. C======================================================================C
  4714.       subroutine xpextr (igrl)
  4715. C                                                                      C
  4716. C     Purpose:  Extrapolate appropriate boundary conditions to the     C
  4717. C               x-plus ( zi-plus ) boundary.                           C
  4718. C                                                                      C
  4719. C======================================================================C
  4720.       include 'UIFlow.com'
  4721.       include 'UIFlow.indx'
  4722. C
  4723.       i = imaxl
  4724. C
  4725.       do 100 n = 1, nsxp
  4726.           jfl  = jfxp(n,igrl)
  4727.       jll  = jlxp(n,igrl)
  4728. C
  4729. C.... Enforce WALL boundary conditions.
  4730. C
  4731.       if ( kbxp(n) .eq. 1 ) then
  4732.        do 10 j   = jfl, jll
  4733.             ij   = ibeg + i + (j-1) * imaxl
  4734.             ij1  = ij - 1
  4735.         p   (ij) = p   (ij1)
  4736.         rho (ij) = rho (ij1)
  4737.         tke (ij) = tke (ij1)
  4738.         tde (ij) = tde (ij1)
  4739.         amut(ij) = amut(ij1)
  4740.         h   (ij) = h   (ij1)
  4741.         t   (ij) = t   (ij1)
  4742.         f   (ij) = f   (ij1)
  4743.         g   (ij) = g   (ij1)
  4744.         yfu (ij) = yfu (ij1)
  4745.         yo2 (ij) = yo2 (ij1)
  4746.         yco2(ij) = yco2(ij1)
  4747.         yh2o(ij) = yh2o(ij1)
  4748.         yn2 (ij) = yn2 (ij1)
  4749. 10     continue
  4750. C
  4751. C.... Enforce INLET boundary condition.
  4752. C
  4753.       elseif ( kbxp(n) .eq. 2 ) then
  4754.        do 20 j = jfl, jll
  4755.            ij  = ibeg + i + (j-1) * imaxl
  4756.            ij1 = ij - 1
  4757.         p(ij)  = p (ij1) + 0.5 * ( p(ij1) - p(ij-2) )
  4758. 20     continue
  4759. C
  4760. C.... Enforce SYMMETRY boundary condition.
  4761. C
  4762.       elseif ( kbxp(n) .eq. 3 ) then
  4763.        do 30 j = jfl, jll
  4764.            ij  = ibeg + i + (j-1) * imaxl
  4765.            ij1 = ij - 1
  4766.         cx   (ij) = 0.0
  4767.         cu   (ij) = 0.0
  4768.         cy   (ij) = cy  (ij1)
  4769.         rho  (ij) = rho (ij1)
  4770.         p    (ij) = p   (ij1)
  4771.         tke  (ij) = tke (ij1)
  4772.         tde  (ij) = tde (ij1)
  4773.         amut (ij) = amut(ij1)
  4774.         h    (ij) = h   (ij1)
  4775.         t    (ij) = t   (ij1)
  4776.         u    (ij) = u   (ij1)
  4777.         v    (ij) = v   (ij1)
  4778.         w    (ij) = w   (ij1)
  4779.         f    (ij) = f   (ij1)
  4780.         g    (ij) = g   (ij1)
  4781.         yfu  (ij) = yfu (ij1)
  4782.         yo2  (ij) = yo2 (ij1)
  4783.         yco2 (ij) = yco2(ij1)
  4784.         yh2o (ij) = yh2o(ij1)
  4785.         yn2  (ij) = yn2 (ij1)
  4786. 30     continue
  4787. C
  4788. C.... Enforce OUTFLOW boundary condition.
  4789. C
  4790.       elseif ( kbxp(n) .eq. 4 ) then
  4791.        do 40 j = jfl, jll
  4792.            ij  = ibeg + i + (j-1) * imaxl
  4793.            ij1 = ij - 1
  4794.         cx   (ij1) = cx (ij-2)
  4795.         cx   (ij)  = cx (ij1)
  4796.         rho  (ij)  = rho(ij1)
  4797.         cu   (ij)  = cx(ij) / rho(ij)
  4798.         cu   (ij1) = cx(ij1) / rho(ij)
  4799.         cy   (ij) = 0.0
  4800.         cv   (ij) = 0.0
  4801.         p    (ij) = p (ij1) + 0.5 * ( p(ij1) - p(ij-2) )
  4802.         tke  (ij) = tke (ij1)
  4803.         tde  (ij) = tde (ij1)
  4804.         amut (ij) = amut(ij1)
  4805.         h    (ij) = h   (ij1)
  4806.         t    (ij) = t   (ij1)
  4807.         u    (ij) = u   (ij1)
  4808.         v    (ij) = 0.0
  4809.         w    (ij) = w   (ij1)
  4810.         f    (ij) = f   (ij1)
  4811.         g    (ij) = g   (ij1)
  4812.         yfu  (ij) = yfu (ij1)
  4813.         yo2  (ij) = yo2 (ij1)
  4814.         yco2 (ij) = yco2(ij1)
  4815.         yh2o (ij) = yh2o(ij1)
  4816.         yn2  (ij) = yn2 (ij1)
  4817. 40     continue
  4818.       endif
  4819. C
  4820. 100   continue
  4821.       return
  4822.       end
  4823. C======================================================================C
  4824.       subroutine ymextr (igrl)
  4825. C                                                                      C
  4826. C     Purpose:  Extrapolate appropriate boundary conditions to the     C
  4827. C               y-minus ( eta-minus ) boundary.                        C
  4828. C                                                                      C
  4829. C======================================================================C
  4830.       include 'UIFlow.com'
  4831.       include 'UIFlow.indx'
  4832. C
  4833.       j = 1
  4834. C
  4835.       do 100 n = 1, nsym
  4836.           ifl  = ifym(n,igrl)
  4837.       ill  = ilym(n,igrl)
  4838. C
  4839. C.... Enforce WALL boundary conditions.
  4840. C
  4841.       if ( kbym(n) .eq. 1 ) then
  4842.        do 10 i   = ifl, ill
  4843.             ij   = ibeg + i + (j-1) * imaxl
  4844.            ij1   = ij + imaxl
  4845.         p   (ij) = p   (ij1)
  4846.         rho (ij) = rho (ij1)
  4847.         tke (ij) = tke (ij1)
  4848.         tde (ij) = tde (ij1)
  4849.         amut(ij) = amut(ij1)
  4850.         h   (ij) = h   (ij1)
  4851.         t   (ij) = t   (ij1)
  4852.         f   (ij) = f   (ij1)
  4853.         g   (ij) = g   (ij1)
  4854.         yfu (ij) = yfu (ij1)
  4855.         yo2 (ij) = yo2 (ij1)
  4856.         yco2(ij) = yco2(ij1)
  4857.         yh2o(ij) = yh2o(ij1)
  4858.         yn2 (ij) = yn2 (ij1)
  4859. 10     continue
  4860. C
  4861. C.... Enforce INLET boundary condition.
  4862. C
  4863.       elseif ( kbym(n) .eq. 2 ) then
  4864.        do 20 i = ifl, ill
  4865.            ij  = ibeg + i + (j-1) * imaxl
  4866.           ij1  = ij + imaxl
  4867.           ij2  = ij1 + imaxl
  4868.         p(ij)  = p (ij1) + 0.5 * ( p(ij1) - p(ij2) )
  4869. 20     continue
  4870. C
  4871. C.... Enforce SYMMETRY boundary condition.
  4872. C
  4873.       elseif ( kbym(n) .eq. 3 ) then
  4874.        do 30 i = ifl, ill
  4875.            ij  = ibeg + i + (j-1) * imaxl
  4876.           ij1  = ij + imaxl
  4877.         cx   (ij) = cx   (ij1)
  4878.         cy   (ij) = 0.0
  4879.         cv   (ij) = 0.0
  4880.         rho  (ij) = rho  (ij1)
  4881.         p    (ij) = p    (ij1)
  4882.         tke  (ij) = tke  (ij1)
  4883.         tde  (ij) = tde  (ij1)
  4884.         amut (ij) = amut (ij1)
  4885.         h    (ij) = h    (ij1)
  4886.         t    (ij) = t    (ij1)
  4887.         u    (ij) = u    (ij1)
  4888.         v    (ij) = 0.0 
  4889.         w    (ij) = 0.0
  4890.         f    (ij) = f    (ij1)
  4891.         g    (ij) = g    (ij1)
  4892.         yfu  (ij) = yfu  (ij1)
  4893.         yo2  (ij) = yo2  (ij1)
  4894.         yh2o (ij) = yh2o (ij1)
  4895.         yco2 (ij) = yco2 (ij1)
  4896.         yn2  (ij) = yn2  (ij1)
  4897. 30     continue
  4898. C
  4899. C.... Enforce OUTFLOW boundary condition.
  4900. C
  4901.       elseif ( kbym(n) .eq. 4 ) then
  4902.        do 40 i = ifl, ill
  4903.            ij  = ibeg + i + (j-1) * imaxl
  4904.           ij1  = ij + imaxl
  4905.           ij2  = ij1 + imaxl
  4906.        cx (ij)   = 0.0
  4907.        cu (ij)   = 0.0
  4908.        rho(ij)   = rho(ij1)
  4909.        cy (ij)   = cy(ij1)
  4910.        cv (ij)   = cy(ij) / rho(ij)
  4911.        p    (ij) = p (ij1) + 0.5 * ( p(ij1) - p(ij2) )
  4912.        tke  (ij) = tke  (ij1)
  4913.        tde  (ij) = tde  (ij1)
  4914.        amut (ij) = amut (ij1)
  4915.        h    (ij) = h    (ij1)
  4916.        t    (ij) = t    (ij1)
  4917.        u    (ij) = u    (ij1)
  4918.        v    (ij) = v    (ij1)
  4919.        w    (ij) = w    (ij1)
  4920.        f    (ij) = f    (ij1)
  4921.        g    (ij) = g    (ij1)
  4922.        yfu  (ij) = yfu  (ij1)
  4923.        yo2  (ij) = yo2  (ij1)
  4924.        yh2o (ij) = yh2o (ij1)
  4925.        yco2 (ij) = yco2 (ij1)
  4926.        yn2  (ij) = yn2  (ij1)
  4927. 40     continue
  4928.       endif
  4929. C
  4930. 100   continue
  4931.       return
  4932.       end
  4933. C======================================================================C
  4934.       subroutine ypextr (igrl)
  4935. C                                                                      C
  4936. C     Purpose:  Extrapolate appropriate boundary conditions to the     C
  4937. C               y-plus ( eta-plus ) boundary.                          C
  4938. C                                                                      C
  4939. C======================================================================C
  4940.       include 'UIFlow.com'
  4941.       include 'UIFlow.indx'
  4942. C
  4943.       j = jmaxl
  4944. C
  4945.       do 100 n = 1, nsyp
  4946.           ifl  = ifyp(n,igrl)
  4947.       ill  = ilyp(n,igrl)
  4948. C
  4949. C.... Enforce WALL boundary conditions.
  4950. C
  4951.       if ( kbyp(n) .eq. 1 ) then
  4952.        do 10 i   = ifl, ill
  4953.             ij   = ibeg + i + (j-1) * imaxl
  4954.            ij1   = ij - imaxl
  4955.         p   (ij) = p   (ij1)
  4956.         rho (ij) = rho (ij1)
  4957.         tke (ij) = tke (ij1)
  4958.         tde (ij) = tde (ij1)
  4959.         amut(ij) = amut(ij1)
  4960.         h   (ij) = h   (ij1)
  4961.         t   (ij) = t   (ij1)
  4962.         f   (ij) = f   (ij1)
  4963.         g   (ij) = g   (ij1)
  4964.         yfu (ij) = yfu (ij1)
  4965.         yo2 (ij) = yo2 (ij1)
  4966.         yco2(ij) = yco2(ij1)
  4967.         yh2o(ij) = yh2o(ij1)
  4968.         yn2 (ij) = yn2 (ij1)
  4969. 10     continue
  4970. C
  4971. C.... Enforce INLET boundary condition.
  4972. C
  4973.       elseif ( kbyp(n) .eq. 2 ) then
  4974.        do 20 i = ifl, ill
  4975.            ij  = ibeg + i + (j-1) * imaxl
  4976.           ij1  = ij - imaxl
  4977.           ij2  = ij1 - imaxl
  4978.         p(ij)  = p (ij1) + 0.5 * ( p(ij1) - p(ij2) )
  4979. 20     continue
  4980. C
  4981. C.... Enforce SYMMETRY boundary condition.
  4982. C
  4983.       elseif ( kbyp(n) .eq. 3 ) then
  4984.        do 30 i = ifl, ill
  4985.            ij  = ibeg + i + (j-1) * imaxl
  4986.           ij1  = ij - imaxl
  4987.         cx   (ij) = cx   (ij1)
  4988.         cy   (ij) = 0.0
  4989.         cv   (ij) = 0.0
  4990.         rho  (ij) = rho  (ij1)
  4991.         p    (ij) = p    (ij1)
  4992.         tke  (ij) = tke  (ij1)
  4993.         tde  (ij) = tde  (ij1)
  4994.         amut (ij) = amut (ij1)
  4995.         h    (ij) = h    (ij1)
  4996.         t    (ij) = t    (ij1)
  4997.         u    (ij) = u    (ij1)
  4998.         v    (ij) = 0.0 
  4999.         w    (ij) = 0.0
  5000.         f    (ij) = f    (ij1)
  5001.         g    (ij) = g    (ij1)
  5002.         yfu  (ij) = yfu  (ij1)
  5003.         yo2  (ij) = yo2  (ij1)
  5004.         yh2o (ij) = yh2o (ij1)
  5005.         yco2 (ij) = yco2 (ij1)
  5006.         yn2  (ij) = yn2  (ij1)
  5007. 30     continue
  5008. C
  5009. C.... Enforce OUTFLOW boundary condition.
  5010. C
  5011.       elseif ( kbyp(n) .eq. 4 ) then
  5012.        do 40 i = ifl, ill
  5013.            ij  = ibeg + i + (j-1) * imaxl
  5014.           ij1  = ij  - imaxl
  5015.           ij2  = ij1 - imaxl
  5016.         cx   (ij) = 0.0
  5017.         cu   (ij) = 0.0
  5018.         rho  (ij) = rho(ij1)
  5019.         cy  (ij1) = cy (ij2)
  5020.         cy   (ij) = cy (ij1)
  5021.         cv  (ij1) = cy (ij1) / rho(ij)
  5022.         cv   (ij) = cv (ij1)
  5023.         p    (ij) = p (ij1) + 0.5 * ( p(ij1) - p(ij2) )
  5024.         tke  (ij) = tke  (ij1)
  5025.         tde  (ij) = tde  (ij1)
  5026.         h    (ij) = h    (ij1)
  5027.         t    (ij) = t    (ij1)
  5028.         u    (ij) = u    (ij1)
  5029.         v    (ij) = v    (ij1)
  5030.         w    (ij) = w    (ij1)
  5031.         f    (ij) = f    (ij1)
  5032.         g    (ij) = g    (ij1)
  5033.         yfu  (ij) = yfu  (ij1)
  5034.         yo2  (ij) = yo2  (ij1)
  5035.         yh2o (ij) = yh2o (ij1)
  5036.         yco2 (ij) = yco2 (ij1)
  5037.         yn2  (ij) = yn2  (ij1)
  5038. 40     continue
  5039.       endif
  5040. C
  5041. 100   continue
  5042.       return
  5043.       end
  5044. C======================================================================C
  5045.       subroutine zero
  5046. C                                                                      C
  5047. C     Purpose:  Set all variables equal to zero.                       C
  5048. C                                                                      C
  5049. C======================================================================C
  5050.       include 'UIFlow.com'
  5051. C
  5052.       do 20 igrl = 1,ngrid
  5053.        nitn(igrl) = 0
  5054.        do 11 nv   =  1, 11
  5055.         error(igrl,nv) = 0.0
  5056. 11     continue
  5057. C
  5058.        if = 1 + nbeg(igrl)
  5059.        il = nbeg(igrl) + imax(igrl) * jmax(igrl)
  5060.        do 10 i = if, il
  5061.         u    (i) = 0.0
  5062.         v    (i) = 0.0
  5063.         w    (i) = 0.0
  5064.         p    (i) = 0.0
  5065.         h    (i) = 0.0
  5066.         t    (i) = 0.0
  5067.         cph  (i) = 0.0
  5068.         pp   (i) = 0.0
  5069.         rho  (i) = 0.0
  5070.         gam  (i) = 0.0
  5071.         tke  (i) = 0.0
  5072.         tde  (i) = 0.0
  5073.         f    (i) = 0.0
  5074.         g    (i) = 0.0
  5075.         yfu  (i) = 0.0
  5076.         yo2  (i) = 0.0
  5077.         yn2  (i) = 0.0
  5078.         yco2 (i) = 0.0
  5079.         yh2o (i) = 0.0
  5080.         x    (i) = 0.0
  5081.         y    (i) = 0.0
  5082.         ajb  (i) = 0.0
  5083.         a11  (i) = 0.0
  5084.         a12  (i) = 0.0
  5085.         a21  (i) = 0.0
  5086.         a22  (i) = 0.0
  5087.         q11  (i) = 0.0
  5088.         q12  (i) = 0.0
  5089.         q21  (i) = 0.0
  5090.         q22  (i) = 0.0
  5091.         fx   (i) = 1.0
  5092.         fy   (i) = 1.0
  5093.         fx1  (i) = 1.0
  5094.         fy1  (i) = 1.0
  5095.         dux  (i) = 0.0
  5096.         duy  (i) = 0.0
  5097.         dvx  (i) = 0.0
  5098.         dvy  (i) = 0.0
  5099.         cu   (i) = 0.0
  5100.         cv   (i) = 0.0
  5101.         cx   (i) = 0.0
  5102.         cy   (i) = 0.0
  5103.         dx   (i) = 0.0
  5104.         dy   (i) = 0.0
  5105.         aw   (i) = 0.0
  5106.         ae   (i) = 0.0
  5107.         as   (i) = 0.0
  5108.         an   (i) = 0.0
  5109.         ap   (i) = 1.0
  5110.         apm  (i) = 1.0
  5111.         app  (i) = 1.0
  5112.         apu  (i) = 0.0
  5113.         apv  (i) = 0.0
  5114.         su   (i) = 0.0
  5115.         sfu  (i) = 0.0
  5116.         sp   (i) = 0.0
  5117.         resu (i) = 0.0
  5118.         resv (i) = 0.0
  5119.         resux(i) = 0.0
  5120.         resvx(i) = 0.0
  5121.         rs   (i) = 0.0
  5122.         dpdx (i) = 0.0
  5123.         dpdy (i) = 0.0
  5124.         amu  (i) = 0.0
  5125.         amut (i) = 0.0
  5126.         prod (i) = 0.0
  5127.         x1   (i) = 0.0
  5128.         x2   (i) = 0.0
  5129.         y1   (i) = 0.0
  5130.         y2   (i) = 0.0
  5131. 10     continue
  5132. 20    continue
  5133. C
  5134.       return
  5135.       end
  5136. C======================================================================C
  5137.       subroutine xwcommon
  5138. C                                                                      C
  5139. C     Purpose:  Set all variables equal to zero.                       C
  5140. C                                                                      C
  5141. C======================================================================C
  5142.       include 'UIFlow.com'
  5143. C
  5144.       open(99, file='TRACE.INPUT' ,access='sequential',status='unknown')
  5145.       write (99, *) 'klam +++'
  5146.       write (99, *) klam,kcomp, kswrl, kpgrid, model
  5147.       write (99, *) 'kfuel +++'
  5148.       write (99, *) kfuel, knorth, kplax, kadj
  5149.       write (99, *) 'ngrid +++'
  5150.       write (99, *) ngrid, ncelx, ncely
  5151. C      
  5152.       write (99, *) '+++ Left Side +++'
  5153.       write (99, *) nsxm
  5154.       do 10 n = 1, nsxm
  5155.         write (99, *) kbxm(n), jfxm, jlxm
  5156.         write (99, *) ubxm(n), vbxm(n), wbxm(n), dvar, txm(n)
  5157.         write (99, *) rhxm(n), fxm(n), gxm(n), tkxm(n), tdxm(n)
  5158.         write (99, *) fuxm(n), co2xm(n), h2oxm(n), o2xm(n), wmxm(n)
  5159. 10    continue
  5160. C      
  5161.       write (99, *) '+++ Right Side +++'
  5162.       write (99, *) nsxp
  5163.       do 20 n = 1,nsxp
  5164.         write (99, *) kbxp (n), jfxp, jlxp
  5165.         write (99, *) ubxp(n), vbxp(n), wbxp(n), dvar, txp(n)
  5166.         write (99, *) rhxp(n), fxp(n), gxp(n), tkxp(n), tdxp(n)
  5167.         write (99, *) fuxp(n), co2xp(n), h2oxp(n), o2xp(n), wmxp(n)
  5168. 20    continue
  5169. C      
  5170.       write (99, *) '+++ Bottom Side +++'
  5171.       write (99, *) nsym
  5172.       do 30 n = 1,nsym
  5173.         write (99, *) kbym (n), ifym, ilym
  5174.         write (99, *) ubym(n), vbym(n), wbym(n), dvar, tym(n)
  5175.         write (99, *) rhym(n), fym(n), gym(n), tkym(n), tdym(n)
  5176.         write (99, *) fuym(n), co2ym(n), h2oym(n), o2ym(n), wmym(n)
  5177. 30    continue
  5178. C      
  5179.       write (99, *) '+++ Top Side +++'
  5180.       write (99, *) nsyp
  5181.       do 40 n = 1,nsyp
  5182.         write (99, *) kbyp(n), ifyp, ilyp
  5183.         write (99, *) ubyp(n), vbyp(n), wbyp(n), dvar, typ(n)
  5184.         write (99, *) rhyp(n), fyp(n), gyp(n), tkyp(n), tdyp(n)
  5185.         write (99, *) fuyp(n), co2yp(n), h2oyp(n), o2yp(n), wmyp(n)
  5186. 40    continue
  5187. C      
  5188.       write (99, *) '+++ Interior +++'
  5189.       write (99, *) ugs, vgs, wgs, rhgs
  5190.       write (99, *) tgs, tkgs, tdgs, fgs, ggs, fugs
  5191.       write (99, *) tfuel, tair
  5192.       write (99, *) (prl(nv), nv=1, 11)
  5193.       write (99, *) (prt(nv), nv=1, 11)
  5194.       write (99, *) (relx(nv), nv=1, 11)
  5195.       write (99, *) (nswp(nv) , nv=1, 11)
  5196.       write (99, *) (iprint(i), i=1 , 12)
  5197.       write (99, *) pref, vscty
  5198.       write (99, *) maxitn, tolr(ngrid)
  5199.       write (99, *) rin
  5200.       write (99, *) nxbaf, nybaf, nobs
  5201. C
  5202.       write (99,*) cd,cdqtr,cdtqtr,ee,ak
  5203.       write (99,*) refu,refv,refw,refc
  5204.       write (99,*) ind1,ind2,ind3,ind4
  5205.       close(UNIT = 99, status='keep')
  5206.       return
  5207.       end
  5208.